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}