Skip to main content

nautilus_model/defi/tick_map/
full_math.rs

1// -------------------------------------------------------------------------------------------------
2//  Copyright (C) 2015-2026 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([u64::from(mm < prod_0), 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([u64::from(remainder > prod_0), 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    #[must_use]
164    pub fn sqrt(x: U256) -> U256 {
165        if x.is_zero() {
166            return U256::ZERO;
167        }
168
169        if x == U256::from(1u128) {
170            return U256::from(1u128);
171        }
172
173        let mut z = x;
174        let mut y = (x + U256::from(1u128)) >> 1;
175
176        while y < z {
177            z = y;
178            y = (x / z + z) >> 1;
179        }
180
181        z
182    }
183
184    /// Truncates a U256 value to u128 by extracting the lower 128 bits.
185    ///
186    /// This matches Solidity's `uint128(value)` cast behavior, which discards
187    /// the upper 128 bits. If the value is larger than `u128::MAX`, the upper
188    /// bits are lost.
189    #[must_use]
190    pub fn truncate_to_u128(value: U256) -> u128 {
191        (value & U256::from(u128::MAX)).to::<u128>()
192    }
193
194    /// Converts an I256 signed integer to U256, mimicking Solidity's `uint256(int256)` cast.
195    ///
196    /// This performs a reinterpret cast, preserving the bit pattern:
197    /// - Positive values: returns the value as-is
198    /// - Negative values: returns the two's complement representation as unsigned
199    #[must_use]
200    pub fn truncate_to_u256(value: I256) -> U256 {
201        value.into_raw()
202    }
203
204    /// Converts a U256 unsigned integer to I256, mimicking Solidity's `int256(uint256)` cast.
205    ///
206    /// This performs a reinterpret cast, preserving the bit pattern.
207    /// Solidity's `SafeCast.toInt256()` checks the value fits in `I256::MAX`, then reinterprets.
208    ///
209    /// # Panics
210    /// Panics if the value exceeds `I256::MAX` (matching Solidity's require check)
211    #[must_use]
212    pub fn truncate_to_i256(value: U256) -> I256 {
213        I256::from_raw(value)
214    }
215}
216
217#[cfg(test)]
218mod tests {
219    use rstest::*;
220
221    use super::*;
222
223    #[rstest]
224    fn test_mul_div_reverts_denominator_zero() {
225        // Test that denominator 0 causes error
226        assert!(FullMath::mul_div(Q128, U256::from(5), U256::ZERO).is_err());
227
228        // Test with numerator overflow and denominator 0
229        assert!(FullMath::mul_div(Q128, Q128, U256::ZERO).is_err());
230    }
231
232    #[rstest]
233    fn test_mul_div_reverts_output_overflow() {
234        // Test output overflow: Q128 * Q128 / 1 would overflow
235        assert!(FullMath::mul_div(Q128, Q128, U256::from(1)).is_err());
236
237        // Test overflow with inputs that would cause prod1 >= denominator
238        // MAX * MAX / 1 would definitely overflow
239        assert!(FullMath::mul_div(U256::MAX, U256::MAX, U256::from(1)).is_err());
240
241        // Test with a smaller denominator that should still cause overflow
242        assert!(FullMath::mul_div(U256::MAX, U256::MAX, U256::from(2)).is_err());
243
244        // Test overflow with all max inputs and denominator = MAX - 1
245        assert!(FullMath::mul_div(U256::MAX, U256::MAX, U256::MAX - U256::from(1)).is_err());
246    }
247
248    #[rstest]
249    fn test_mul_div_all_max_inputs() {
250        // MAX * MAX / MAX = MAX
251        let result = FullMath::mul_div(U256::MAX, U256::MAX, U256::MAX).unwrap();
252        assert_eq!(result, U256::MAX);
253    }
254
255    #[rstest]
256    fn test_mul_div_accurate_without_phantom_overflow() {
257        // Calculate Q128 * 0.5 / 1.5 = Q128 / 3
258        let numerator_b = Q128 * U256::from(50) / U256::from(100); // 0.5
259        let denominator = Q128 * U256::from(150) / U256::from(100); // 1.5
260        let expected = Q128 / U256::from(3);
261
262        let result = FullMath::mul_div(Q128, numerator_b, denominator).unwrap();
263        assert_eq!(result, expected);
264    }
265
266    #[rstest]
267    fn test_mul_div_accurate_with_phantom_overflow() {
268        // Calculate Q128 * 35 * Q128 / (8 * Q128) = 35/8 * Q128 = 4.375 * Q128
269        let numerator_b = U256::from(35) * Q128;
270        let denominator = U256::from(8) * Q128;
271        let expected = U256::from(4375) * Q128 / U256::from(1000);
272
273        let result = FullMath::mul_div(Q128, numerator_b, denominator).unwrap();
274        assert_eq!(result, expected);
275    }
276
277    #[rstest]
278    fn test_mul_div_accurate_with_phantom_overflow_repeating_decimal() {
279        // Calculate Q128 * 1000 * Q128 / (3000 * Q128) = 1/3 * Q128
280        let numerator_b = U256::from(1000) * Q128;
281        let denominator = U256::from(3000) * Q128;
282        let expected = Q128 / U256::from(3);
283
284        let result = FullMath::mul_div(Q128, numerator_b, denominator).unwrap();
285        assert_eq!(result, expected);
286    }
287
288    #[rstest]
289    fn test_mul_div_basic_cases() {
290        // Simple case: 100 * 200 / 50 = 400
291        assert_eq!(
292            FullMath::mul_div(U256::from(100), U256::from(200), U256::from(50)).unwrap(),
293            U256::from(400)
294        );
295
296        // Test with 1: a * 1 / b = a / b
297        assert_eq!(
298            FullMath::mul_div(U256::from(1000), U256::from(1), U256::from(4)).unwrap(),
299            U256::from(250)
300        );
301
302        // Test division that results in 0 due to floor
303        assert_eq!(
304            FullMath::mul_div(U256::from(1), U256::from(1), U256::from(3)).unwrap(),
305            U256::ZERO
306        );
307    }
308
309    // mul_div_rounding_up tests
310    #[rstest]
311    fn test_mul_div_rounding_up_reverts_denominator_zero() {
312        // Test that denominator 0 causes error
313        assert!(FullMath::mul_div_rounding_up(Q128, U256::from(5), U256::ZERO).is_err());
314
315        // Test with numerator overflow and denominator 0
316        assert!(FullMath::mul_div_rounding_up(Q128, Q128, U256::ZERO).is_err());
317    }
318
319    #[rstest]
320    fn test_mul_div_rounding_up_reverts_output_overflow() {
321        // Test output overflow: Q128 * Q128 / 1 would overflow
322        assert!(FullMath::mul_div_rounding_up(Q128, Q128, U256::from(1)).is_err());
323
324        // Test overflow with all max inputs minus 1 - this should pass since MAX/MAX-1 = ~1
325        // but since there's a remainder, rounding up would still fit in U256
326        // Let's test a case that actually overflows after rounding
327        assert!(FullMath::mul_div_rounding_up(U256::MAX, U256::MAX, U256::from(2)).is_err());
328
329        // Test overflow with all max inputs and denominator = MAX - 1
330        assert!(
331            FullMath::mul_div_rounding_up(U256::MAX, U256::MAX, U256::MAX - U256::from(1)).is_err()
332        );
333    }
334
335    #[rstest]
336    fn test_mul_div_rounding_up_reverts_overflow_after_rounding_case_1() {
337        // Edge case discovered through fuzzing: mul_div succeeds but result is MAX with remainder
338        // so rounding up would overflow
339        let a = U256::from_str_radix("535006138814359", 10).unwrap();
340        let b = U256::from_str_radix(
341            "432862656469423142931042426214547535783388063929571229938474969",
342            10,
343        )
344        .unwrap();
345        let denominator = U256::from(2);
346
347        assert!(FullMath::mul_div_rounding_up(a, b, denominator).is_err());
348    }
349
350    #[rstest]
351    fn test_mul_div_rounding_up_reverts_overflow_after_rounding_case_2() {
352        // Another edge case discovered through fuzzing: tests boundary condition where
353        // mul_div returns MAX-1 but with remainder, so rounding up would cause overflow
354        let a = U256::from_str_radix(
355            "115792089237316195423570985008687907853269984659341747863450311749907997002549",
356            10,
357        )
358        .unwrap();
359        let b = U256::from_str_radix(
360            "115792089237316195423570985008687907853269984659341747863450311749907997002550",
361            10,
362        )
363        .unwrap();
364        let denominator = U256::from_str_radix(
365            "115792089237316195423570985008687907853269984653042931687443039491902864365164",
366            10,
367        )
368        .unwrap();
369
370        assert!(FullMath::mul_div_rounding_up(a, b, denominator).is_err());
371    }
372
373    #[rstest]
374    fn test_mul_div_rounding_up_all_max_inputs() {
375        // MAX * MAX / MAX = MAX (no rounding needed)
376        let result = FullMath::mul_div_rounding_up(U256::MAX, U256::MAX, U256::MAX).unwrap();
377        assert_eq!(result, U256::MAX);
378    }
379
380    #[rstest]
381    fn test_mul_div_rounding_up_accurate_without_phantom_overflow() {
382        // Calculate Q128 * 0.5 / 1.5 = Q128 / 3, but with rounding up
383        let numerator_b = Q128 * U256::from(50) / U256::from(100); // 0.5
384        let denominator = Q128 * U256::from(150) / U256::from(100); // 1.5
385        let expected = Q128 / U256::from(3) + U256::from(1); // Rounded up
386
387        let result = FullMath::mul_div_rounding_up(Q128, numerator_b, denominator).unwrap();
388        assert_eq!(result, expected);
389    }
390
391    #[rstest]
392    fn test_mul_div_rounding_up_accurate_with_phantom_overflow() {
393        // Calculate Q128 * 35 * Q128 / (8 * Q128) = 35/8 * Q128 = 4.375 * Q128
394        // This should be exact (no remainder), so no rounding up needed
395        let numerator_b = U256::from(35) * Q128;
396        let denominator = U256::from(8) * Q128;
397        let expected = U256::from(4375) * Q128 / U256::from(1000);
398
399        let result = FullMath::mul_div_rounding_up(Q128, numerator_b, denominator).unwrap();
400        assert_eq!(result, expected);
401    }
402
403    #[rstest]
404    fn test_mul_div_rounding_up_accurate_with_phantom_overflow_repeating_decimal() {
405        // Calculate Q128 * 1000 * Q128 / (3000 * Q128) = 1/3 * Q128, with rounding up
406        let numerator_b = U256::from(1000) * Q128;
407        let denominator = U256::from(3000) * Q128;
408        let expected = Q128 / U256::from(3) + U256::from(1); // Rounded up due to remainder
409
410        let result = FullMath::mul_div_rounding_up(Q128, numerator_b, denominator).unwrap();
411        assert_eq!(result, expected);
412    }
413
414    #[rstest]
415    fn test_mul_div_rounding_up_basic_cases() {
416        // Test exact division (no rounding needed)
417        assert_eq!(
418            FullMath::mul_div_rounding_up(U256::from(100), U256::from(200), U256::from(50))
419                .unwrap(),
420            U256::from(400)
421        );
422
423        // Test division with remainder (rounding up needed)
424        assert_eq!(
425            FullMath::mul_div_rounding_up(U256::from(1), U256::from(1), U256::from(3)).unwrap(),
426            U256::from(1) // 0 rounded up to 1
427        );
428
429        // Test another rounding case: 7 * 3 / 4 = 21 / 4 = 5.25 -> 6
430        assert_eq!(
431            FullMath::mul_div_rounding_up(U256::from(7), U256::from(3), U256::from(4)).unwrap(),
432            U256::from(6)
433        );
434
435        // Test case with zero result and zero remainder
436        assert_eq!(
437            FullMath::mul_div_rounding_up(U256::ZERO, U256::from(100), U256::from(3)).unwrap(),
438            U256::ZERO
439        );
440    }
441
442    #[rstest]
443    fn test_mul_div_rounding_up_overflow_at_max() {
444        // Test that rounding up when result is already MAX causes overflow
445        // We need a case where mul_div returns MAX but there's a remainder
446        // This is tricky to construct, so we test the boundary condition
447        assert!(FullMath::mul_div_rounding_up(U256::MAX, U256::from(2), U256::from(2)).is_ok());
448
449        // This should succeed: MAX * 1 / 1 = MAX (no remainder)
450        assert_eq!(
451            FullMath::mul_div_rounding_up(U256::MAX, U256::from(1), U256::from(1)).unwrap(),
452            U256::MAX
453        );
454    }
455
456    #[rstest]
457    fn test_truncate_to_u128_preserves_small_values() {
458        // Small value (fits in u128) should be preserved exactly
459        let value = U256::from(12345u128);
460        assert_eq!(FullMath::truncate_to_u128(value), 12345u128);
461
462        // u128::MAX should be preserved
463        let max_value = U256::from(u128::MAX);
464        assert_eq!(FullMath::truncate_to_u128(max_value), u128::MAX);
465    }
466
467    #[rstest]
468    fn test_truncate_to_u128_discards_upper_bits() {
469        // Value = u128::MAX + 1 (sets bit 128)
470        // Lower 128 bits = 0, so result should be 0
471        let value = U256::from(u128::MAX) + U256::from(1);
472        assert_eq!(FullMath::truncate_to_u128(value), 0);
473
474        // Value with both high and low bits set:
475        // High 128 bits = 0xFFFF...FFFF, Low 128 bits = 0x1234
476        let value = (U256::from(u128::MAX) << 128) | U256::from(0x1234u128);
477        assert_eq!(FullMath::truncate_to_u128(value), 0x1234u128);
478    }
479}