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}