nautilus_model/defi/tick_map/
full_math.rs

1// -------------------------------------------------------------------------------------------------
2//  Copyright (C) 2015-2025 Nautech Systems Pty Ltd. All rights reserved.
3//  https://nautechsystems.io
4//
5//  Licensed under the GNU Lesser General Public License Version 3.0 (the "License");
6//  You may not use this file except in compliance with the License.
7//  You may obtain a copy of the License at https://www.gnu.org/licenses/lgpl-3.0.en.html
8//
9//  Unless required by applicable law or agreed to in writing, software
10//  distributed under the License is distributed on an "AS IS" BASIS,
11//  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12//  See the License for the specific language governing permissions and
13//  limitations under the License.
14// -------------------------------------------------------------------------------------------------
15
16use alloy_primitives::{I256, U160, U256};
17
18pub const Q128: U256 = U256::from_limbs([0, 0, 1, 0]);
19pub const Q96_U160: U160 = U160::from_limbs([0, 1 << 32, 0]);
20
21/// Contains 512-bit math functions for Uniswap V3 style calculations
22/// Handles "phantom overflow" - allows multiplication and division where
23/// intermediate values overflow 256 bits
24#[derive(Debug)]
25pub struct FullMath;
26
27impl FullMath {
28    /// Calculates floor(a×b÷denominator) with full precision
29    ///
30    /// Follows the Solidity implementation from Uniswap V3's FullMath library:
31    /// <https://github.com/Uniswap/v3-core/blob/main/contracts/libraries/FullMath.sol>
32    ///
33    /// # Errors
34    ///
35    /// Returns error if `denominator` is zero or the result would overflow 256 bits.
36    pub fn mul_div(a: U256, b: U256, mut denominator: U256) -> anyhow::Result<U256> {
37        // 512-bit multiply [prod1 prod0] = a * b
38        // Compute the product mod 2**256 and mod 2**256 - 1
39        // then use the Chinese Remainder Theorem to reconstruct
40        // the 512 bit result. The result is stored in two 256
41        // variables such that product = prod1 * 2**256 + prod0
42        let mm = a.mul_mod(b, U256::MAX);
43
44        // Least significant 256 bits of the product
45        let mut prod_0 = a * b;
46        let mut prod_1 = mm - prod_0 - U256::from_limbs([(mm < prod_0) as u64, 0, 0, 0]);
47
48        // Make sure the result is less than 2**256.
49        // Also prevents denominator == 0
50        if denominator <= prod_1 {
51            anyhow::bail!("Result would overflow 256 bits");
52        }
53
54        ///////////////////////////////////////////////
55        // 512 by 256 division.
56        ///////////////////////////////////////////////
57
58        // Make division exact by subtracting the remainder from [prod1 prod0]
59        // Compute remainder using mul_mod
60        let remainder = a.mul_mod(b, denominator);
61
62        // Subtract 256 bit number from 512 bit number
63        prod_1 -= U256::from_limbs([(remainder > prod_0) as u64, 0, 0, 0]);
64        prod_0 -= remainder;
65
66        // Factor powers of two out of denominator
67        // Compute largest power of two divisor of denominator.
68        // Always >= 1.
69        let mut twos = (-denominator) & denominator;
70
71        // Divide denominator by power of two
72        denominator /= twos;
73
74        // Divide [prod1 prod0] by the factors of two
75        prod_0 /= twos;
76
77        // Shift in bits from prod1 into prod0. For this we need
78        // to flip `twos` such that it is 2**256 / twos.
79        // If twos is zero, then it becomes one
80        twos = (-twos) / twos + U256::from(1);
81
82        prod_0 |= prod_1 * twos;
83
84        // Invert denominator mod 2**256
85        // Now that denominator is an odd number, it has an inverse
86        // modulo 2**256 such that denominator * inv = 1 mod 2**256.
87        // Compute the inverse by starting with a seed that is correct
88        // correct for four bits. That is, denominator * inv = 1 mod 2**4
89        let mut inv = (U256::from(3) * denominator) ^ U256::from(2);
90
91        // Now use Newton-Raphson iteration to improve the precision.
92        // Thanks to Hensel's lifting lemma, this also works in modular
93        // arithmetic, doubling the correct bits in each step.
94        inv *= U256::from(2) - denominator * inv; // inverse mod 2**8
95
96        inv *= U256::from(2) - denominator * inv; // inverse mod 2**16
97
98        inv *= U256::from(2) - denominator * inv; // inverse mod 2**32
99
100        inv *= U256::from(2) - denominator * inv; // inverse mod 2**64
101
102        inv *= U256::from(2) - denominator * inv; // inverse mod 2**128
103
104        inv *= U256::from(2) - denominator * inv; // inverse mod 2**256
105
106        // Because the division is now exact we can divide by multiplying
107        // with the modular inverse of denominator. This will give us the
108        // correct result modulo 2**256. Since the preconditions guarantee
109        // that the outcome is less than 2**256, this is the final result.
110        // We don't need to compute the high bits of the result and prod1
111        // is no longer required.
112        let result = prod_0 * inv;
113
114        Ok(result)
115    }
116
117    /// Calculates ceil(a×b÷denominator) with full precision
118    /// Returns `Ok` with the rounded result or an error when rounding cannot be performed safely.
119    ///
120    /// # Errors
121    ///
122    /// Returns error if `denominator` is zero or the rounded result would overflow `U256`.
123    pub fn mul_div_rounding_up(a: U256, b: U256, denominator: U256) -> anyhow::Result<U256> {
124        let result = Self::mul_div(a, b, denominator)?;
125
126        // Check if there's a remainder
127        if a.mul_mod(b, denominator).is_zero() {
128            Ok(result)
129        } else if result == U256::MAX {
130            anyhow::bail!("Result would overflow 256 bits")
131        } else {
132            Ok(result + U256::from(1))
133        }
134    }
135
136    /// Calculates ceil(a÷b) with proper rounding up
137    /// Equivalent to Solidity's divRoundingUp function
138    ///
139    /// # Errors
140    ///
141    /// Returns error if `b` is zero or if the rounded quotient would overflow `U256`.
142    pub fn div_rounding_up(a: U256, b: U256) -> anyhow::Result<U256> {
143        if b.is_zero() {
144            anyhow::bail!("Cannot divide by zero");
145        }
146
147        let quotient = a / b;
148        let remainder = a % b;
149
150        // Add 1 if there's a remainder (equivalent to gt(mod(x, y), 0) in assembly)
151        if remainder > U256::ZERO {
152            // Check for overflow before incrementing
153            if quotient == U256::MAX {
154                anyhow::bail!("Result would overflow 256 bits");
155            }
156            Ok(quotient + U256::from(1))
157        } else {
158            Ok(quotient)
159        }
160    }
161
162    /// Computes the integer square root of a 256-bit unsigned integer using the Babylonian method
163    pub fn sqrt(x: U256) -> U256 {
164        if x.is_zero() {
165            return U256::ZERO;
166        }
167        if x == U256::from(1u128) {
168            return U256::from(1u128);
169        }
170
171        let mut z = x;
172        let mut y = (x + U256::from(1u128)) >> 1;
173
174        while y < z {
175            z = y;
176            y = (x / z + z) >> 1;
177        }
178
179        z
180    }
181
182    /// Truncates a U256 value to u128 by extracting the lower 128 bits.
183    ///
184    /// This matches Solidity's `uint128(value)` cast behavior, which discards
185    /// the upper 128 bits. If the value is larger than `u128::MAX`, the upper
186    /// bits are lost.
187    #[must_use]
188    pub fn truncate_to_u128(value: U256) -> u128 {
189        (value & U256::from(u128::MAX)).to::<u128>()
190    }
191
192    /// Converts an I256 signed integer to U256, mimicking Solidity's `uint256(int256)` cast.
193    ///
194    /// This performs a reinterpret cast, preserving the bit pattern:
195    /// - Positive values: returns the value as-is
196    /// - Negative values: returns the two's complement representation as unsigned
197    #[must_use]
198    pub fn truncate_to_u256(value: I256) -> U256 {
199        value.into_raw()
200    }
201
202    /// Converts a U256 unsigned integer to I256, mimicking Solidity's `int256(uint256)` cast.
203    ///
204    /// This performs a reinterpret cast, preserving the bit pattern.
205    /// Solidity's SafeCast.toInt256() checks the value fits in I256::MAX, then reinterprets.
206    ///
207    /// # Panics
208    /// Panics if the value exceeds I256::MAX (matching Solidity's require check)
209    #[must_use]
210    pub fn truncate_to_i256(value: U256) -> I256 {
211        I256::from_raw(value)
212    }
213}
214
215#[cfg(test)]
216mod tests {
217    use rstest::*;
218
219    use super::*;
220
221    #[rstest]
222    fn test_mul_div_reverts_denominator_zero() {
223        // Test that denominator 0 causes error
224        assert!(FullMath::mul_div(Q128, U256::from(5), U256::ZERO).is_err());
225
226        // Test with numerator overflow and denominator 0
227        assert!(FullMath::mul_div(Q128, Q128, U256::ZERO).is_err());
228    }
229
230    #[rstest]
231    fn test_mul_div_reverts_output_overflow() {
232        // Test output overflow: Q128 * Q128 / 1 would overflow
233        assert!(FullMath::mul_div(Q128, Q128, U256::from(1)).is_err());
234
235        // Test overflow with inputs that would cause prod1 >= denominator
236        // MAX * MAX / 1 would definitely overflow
237        assert!(FullMath::mul_div(U256::MAX, U256::MAX, U256::from(1)).is_err());
238
239        // Test with a smaller denominator that should still cause overflow
240        assert!(FullMath::mul_div(U256::MAX, U256::MAX, U256::from(2)).is_err());
241
242        // Test overflow with all max inputs and denominator = MAX - 1
243        assert!(FullMath::mul_div(U256::MAX, U256::MAX, U256::MAX - U256::from(1)).is_err());
244    }
245
246    #[rstest]
247    fn test_mul_div_all_max_inputs() {
248        // MAX * MAX / MAX = MAX
249        let result = FullMath::mul_div(U256::MAX, U256::MAX, U256::MAX).unwrap();
250        assert_eq!(result, U256::MAX);
251    }
252
253    #[rstest]
254    fn test_mul_div_accurate_without_phantom_overflow() {
255        // Calculate Q128 * 0.5 / 1.5 = Q128 / 3
256        let numerator_b = Q128 * U256::from(50) / U256::from(100); // 0.5
257        let denominator = Q128 * U256::from(150) / U256::from(100); // 1.5
258        let expected = Q128 / U256::from(3);
259
260        let result = FullMath::mul_div(Q128, numerator_b, denominator).unwrap();
261        assert_eq!(result, expected);
262    }
263
264    #[rstest]
265    fn test_mul_div_accurate_with_phantom_overflow() {
266        // Calculate Q128 * 35 * Q128 / (8 * Q128) = 35/8 * Q128 = 4.375 * Q128
267        let numerator_b = U256::from(35) * Q128;
268        let denominator = U256::from(8) * Q128;
269        let expected = U256::from(4375) * Q128 / U256::from(1000);
270
271        let result = FullMath::mul_div(Q128, numerator_b, denominator).unwrap();
272        assert_eq!(result, expected);
273    }
274
275    #[rstest]
276    fn test_mul_div_accurate_with_phantom_overflow_repeating_decimal() {
277        // Calculate Q128 * 1000 * Q128 / (3000 * Q128) = 1/3 * Q128
278        let numerator_b = U256::from(1000) * Q128;
279        let denominator = U256::from(3000) * Q128;
280        let expected = Q128 / U256::from(3);
281
282        let result = FullMath::mul_div(Q128, numerator_b, denominator).unwrap();
283        assert_eq!(result, expected);
284    }
285
286    #[rstest]
287    fn test_mul_div_basic_cases() {
288        // Simple case: 100 * 200 / 50 = 400
289        assert_eq!(
290            FullMath::mul_div(U256::from(100), U256::from(200), U256::from(50)).unwrap(),
291            U256::from(400)
292        );
293
294        // Test with 1: a * 1 / b = a / b
295        assert_eq!(
296            FullMath::mul_div(U256::from(1000), U256::from(1), U256::from(4)).unwrap(),
297            U256::from(250)
298        );
299
300        // Test division that results in 0 due to floor
301        assert_eq!(
302            FullMath::mul_div(U256::from(1), U256::from(1), U256::from(3)).unwrap(),
303            U256::ZERO
304        );
305    }
306
307    // mul_div_rounding_up tests
308    #[rstest]
309    fn test_mul_div_rounding_up_reverts_denominator_zero() {
310        // Test that denominator 0 causes error
311        assert!(FullMath::mul_div_rounding_up(Q128, U256::from(5), U256::ZERO).is_err());
312
313        // Test with numerator overflow and denominator 0
314        assert!(FullMath::mul_div_rounding_up(Q128, Q128, U256::ZERO).is_err());
315    }
316
317    #[rstest]
318    fn test_mul_div_rounding_up_reverts_output_overflow() {
319        // Test output overflow: Q128 * Q128 / 1 would overflow
320        assert!(FullMath::mul_div_rounding_up(Q128, Q128, U256::from(1)).is_err());
321
322        // Test overflow with all max inputs minus 1 - this should pass since MAX/MAX-1 = ~1
323        // but since there's a remainder, rounding up would still fit in U256
324        // Let's test a case that actually overflows after rounding
325        assert!(FullMath::mul_div_rounding_up(U256::MAX, U256::MAX, U256::from(2)).is_err());
326
327        // Test overflow with all max inputs and denominator = MAX - 1
328        assert!(
329            FullMath::mul_div_rounding_up(U256::MAX, U256::MAX, U256::MAX - U256::from(1)).is_err()
330        );
331    }
332
333    #[rstest]
334    fn test_mul_div_rounding_up_reverts_overflow_after_rounding_case_1() {
335        // Edge case discovered through fuzzing: mul_div succeeds but result is MAX with remainder
336        // so rounding up would overflow
337        let a = U256::from_str_radix("535006138814359", 10).unwrap();
338        let b = U256::from_str_radix(
339            "432862656469423142931042426214547535783388063929571229938474969",
340            10,
341        )
342        .unwrap();
343        let denominator = U256::from(2);
344
345        assert!(FullMath::mul_div_rounding_up(a, b, denominator).is_err());
346    }
347
348    #[rstest]
349    fn test_mul_div_rounding_up_reverts_overflow_after_rounding_case_2() {
350        // Another edge case discovered through fuzzing: tests boundary condition where
351        // mul_div returns MAX-1 but with remainder, so rounding up would cause overflow
352        let a = U256::from_str_radix(
353            "115792089237316195423570985008687907853269984659341747863450311749907997002549",
354            10,
355        )
356        .unwrap();
357        let b = U256::from_str_radix(
358            "115792089237316195423570985008687907853269984659341747863450311749907997002550",
359            10,
360        )
361        .unwrap();
362        let denominator = U256::from_str_radix(
363            "115792089237316195423570985008687907853269984653042931687443039491902864365164",
364            10,
365        )
366        .unwrap();
367
368        assert!(FullMath::mul_div_rounding_up(a, b, denominator).is_err());
369    }
370
371    #[rstest]
372    fn test_mul_div_rounding_up_all_max_inputs() {
373        // MAX * MAX / MAX = MAX (no rounding needed)
374        let result = FullMath::mul_div_rounding_up(U256::MAX, U256::MAX, U256::MAX).unwrap();
375        assert_eq!(result, U256::MAX);
376    }
377
378    #[rstest]
379    fn test_mul_div_rounding_up_accurate_without_phantom_overflow() {
380        // Calculate Q128 * 0.5 / 1.5 = Q128 / 3, but with rounding up
381        let numerator_b = Q128 * U256::from(50) / U256::from(100); // 0.5
382        let denominator = Q128 * U256::from(150) / U256::from(100); // 1.5
383        let expected = Q128 / U256::from(3) + U256::from(1); // Rounded up
384
385        let result = FullMath::mul_div_rounding_up(Q128, numerator_b, denominator).unwrap();
386        assert_eq!(result, expected);
387    }
388
389    #[rstest]
390    fn test_mul_div_rounding_up_accurate_with_phantom_overflow() {
391        // Calculate Q128 * 35 * Q128 / (8 * Q128) = 35/8 * Q128 = 4.375 * Q128
392        // This should be exact (no remainder), so no rounding up needed
393        let numerator_b = U256::from(35) * Q128;
394        let denominator = U256::from(8) * Q128;
395        let expected = U256::from(4375) * Q128 / U256::from(1000);
396
397        let result = FullMath::mul_div_rounding_up(Q128, numerator_b, denominator).unwrap();
398        assert_eq!(result, expected);
399    }
400
401    #[rstest]
402    fn test_mul_div_rounding_up_accurate_with_phantom_overflow_repeating_decimal() {
403        // Calculate Q128 * 1000 * Q128 / (3000 * Q128) = 1/3 * Q128, with rounding up
404        let numerator_b = U256::from(1000) * Q128;
405        let denominator = U256::from(3000) * Q128;
406        let expected = Q128 / U256::from(3) + U256::from(1); // Rounded up due to remainder
407
408        let result = FullMath::mul_div_rounding_up(Q128, numerator_b, denominator).unwrap();
409        assert_eq!(result, expected);
410    }
411
412    #[rstest]
413    fn test_mul_div_rounding_up_basic_cases() {
414        // Test exact division (no rounding needed)
415        assert_eq!(
416            FullMath::mul_div_rounding_up(U256::from(100), U256::from(200), U256::from(50))
417                .unwrap(),
418            U256::from(400)
419        );
420
421        // Test division with remainder (rounding up needed)
422        assert_eq!(
423            FullMath::mul_div_rounding_up(U256::from(1), U256::from(1), U256::from(3)).unwrap(),
424            U256::from(1) // 0 rounded up to 1
425        );
426
427        // Test another rounding case: 7 * 3 / 4 = 21 / 4 = 5.25 -> 6
428        assert_eq!(
429            FullMath::mul_div_rounding_up(U256::from(7), U256::from(3), U256::from(4)).unwrap(),
430            U256::from(6)
431        );
432
433        // Test case with zero result and zero remainder
434        assert_eq!(
435            FullMath::mul_div_rounding_up(U256::ZERO, U256::from(100), U256::from(3)).unwrap(),
436            U256::ZERO
437        );
438    }
439
440    #[rstest]
441    fn test_mul_div_rounding_up_overflow_at_max() {
442        // Test that rounding up when result is already MAX causes overflow
443        // We need a case where mul_div returns MAX but there's a remainder
444        // This is tricky to construct, so we test the boundary condition
445        assert!(FullMath::mul_div_rounding_up(U256::MAX, U256::from(2), U256::from(2)).is_ok());
446
447        // This should succeed: MAX * 1 / 1 = MAX (no remainder)
448        assert_eq!(
449            FullMath::mul_div_rounding_up(U256::MAX, U256::from(1), U256::from(1)).unwrap(),
450            U256::MAX
451        );
452    }
453
454    #[rstest]
455    fn test_truncate_to_u128_preserves_small_values() {
456        // Small value (fits in u128) should be preserved exactly
457        let value = U256::from(12345u128);
458        assert_eq!(FullMath::truncate_to_u128(value), 12345u128);
459
460        // u128::MAX should be preserved
461        let max_value = U256::from(u128::MAX);
462        assert_eq!(FullMath::truncate_to_u128(max_value), u128::MAX);
463    }
464
465    #[rstest]
466    fn test_truncate_to_u128_discards_upper_bits() {
467        // Value = u128::MAX + 1 (sets bit 128)
468        // Lower 128 bits = 0, so result should be 0
469        let value = U256::from(u128::MAX) + U256::from(1);
470        assert_eq!(FullMath::truncate_to_u128(value), 0);
471
472        // Value with both high and low bits set:
473        // High 128 bits = 0xFFFF...FFFF, Low 128 bits = 0x1234
474        let value = (U256::from(u128::MAX) << 128) | U256::from(0x1234u128);
475        assert_eq!(FullMath::truncate_to_u128(value), 0x1234u128);
476    }
477}