Deep Diving into PRBMath
There are severals reasons why doing math is hard in Solidity, with the chief one being the lack of support for fixed-point types. As per the documentation, fixed-point numbers can be declared, but not assigned to or from. Developers are expected to figure out on their own how to implement advanced math functions like exponentials and logarithms.
PRBMath is a smart contract library that purports to solve the aforementioned difficulties and make math in Solidity easy and fun.
Overview
PRBMath comes in four flavors:
- PRBMathSD59x18
- PRBMathSD59x18Typed
- PRBMathUD60x18
- PRBMathUD60x18Typed
The first two work with signed 59.18-decimal fixed-point numbers, while the latter two work with unsigned 60.18-decimal fixed-point numbers. The first and the third operate on basic integers, while the second and the fourth operate on structs. There is no qualitative difference between the non-typed and the typed versions - you use one over the other when you want to demarcate between basic integers and fixed-point numbers in your contracts.
The values the numbers can take are bound by the minimum and the maximum values permitted by the Solidity types int256 and uint256, while the name of the number formats stems from the fact that there can be up to 59/60 digits in the integer part and up to 18 decimals in the fractional part.
These are the functions that you are getting out of the box:
- Absolute
- Arithmetic and geometric average
- Exponentials (binary and natural)
- Floor and ceil
- Fractional
- Inverse
- Logarithms (binary, common and natural)
- Powers (fractional number and basic integers as exponents)
- Multiplication and division
- Square root
In addition, there are getters for mathematical constants:
- Euler’s number
- Pi
- Scale (1e18, which is 1 in fixed-point representation)
Finally, there are functions to convert from and to the fixed-point representation.
How to Use
The README on GitHub offers comprehensive instructions for how to use PRBMath, but here’s a sample:
// SPDX-License-Identifier: UNLICENSED
pragma solidity >=0.8.0;
import "prb-math/contracts/PRBMathSD59x18.sol";
contract SignedConsumer {
using PRBMathSD59x18 for int256;
function foo(int256 x, int256 y) external pure {
int256 result;
result = x.abs();
result = x.avg(y);
result = x.ceil();
result = x.div(y);
result = x.e();
result = x.exp();
result = x.exp2();
result = x.floor();
result = x.frac();
result = x.gm(y);
result = x.inv();
result = x.ln();
result = x.log10();
result = x.log2();
result = x.pi();
result = x.mul(y);
result = x.pow(2);
result = x.scale();
result = x.sqrt();
}
}
To make the code compile, you will need to use either yarn or npm to install the prb-math package.
Target Audience
There’s not much I can predict in terms of what the library is going to be used for, specifically. If anything, we will make good use of it in the next release of our Hifi protocol.
What I can think of though are the following use cases:
- Compounding interest rates.
- TWAP oracles.
- Pricing derivatives like options.
- Standalone AMMs.
- Custom price curves for Uniswap v3 and Balancer v2.
Parts and Parcels
PRBMath is the result of a month-long tinkering with math puzzles and computer science algorithms. My goal was to write a math library that is at the same time practical, intuitive and efficient.
To obtain the best possible gas results, I resorted to assembly code and unchecked arithmetic. I was surprised to notice that in many cases it’s cheaper to check for overflow manually than let the Solidity compiler do it via checked arithmetic.
Let’s delve into each function. Note that I will use the non-typed SD59x18 flavour throughout this document, but you can substite it for UD60x18 in all cases and the meaning remains the same. Further note that I will make frequent use of the “phantom overflow” notion, which refers to certain mathematical operations whereby the intermediary result overflows but the final result doesn’t.
Abs
Signature
function abs(int256 x) internal pure returns (int256 result);
Explanation
This is a simple function. The first thing I do is check that x is not min sd59x18, since in this case the absolute value couldn’t be represented in Solidity (read about two’s complement). Then I use the unary operator -
to flip the value if x is negative or simply return the value if x is positive.
Avg
Signature
function avg(int256 x, int256 y) internal pure returns (int256 result);
Explanation
The naive way of calculating the average of two numbers is to add them together, then divide by two. The issue with this approach is that it’s susceptible to phantom overflow.
To go around this, I divide each number by two and then add the result. But I must not forget that Solidity truncates the result when performing division. Thus, if both inputs are odd, I have to add one to the result to account for the doubly truncated halves.
Bonus: the chronicle of ideas that I iterated through before arriving at the fastest solution.
Ceil
Signature
function ceil(int256 x) internal pure returns (int256 result);
Explanation
I confess that I didn’t expect this function to be as subtle to implement as it was. This fact that ethers handles the modulo function differently to Solidity only made matters worse.
The first thing I do is check that x is less than or equal to the value dubbed “max whole sd59x18”. This is basically max sd59x18 but with all 18 decimals set to zero. This requirement is in place because the ceil couldn’t otherwise be represented in Solidity.
Then, the idea is to calculate the remainder of x divided by the scale number (1e18), and:
- If the remainder is zero, return x.
- If the remainder is non-zero, I set the result to the difference between x and the remainder calculated above. And finally, if x is positive, I also add the scale number to the result.
Div
Signature
function div(int256 x, int256 y) internal pure returns (int256 result);
Explanation
Dividing two fixed-point numbers is done by multiplying the first fixed-point number by the scale number 1e18, and then dividing the product by second fixed-point number.
Naturally, this creates a risk of phantom overflow. I prevent this by using Remco Bloemen’s mathematical mulDiv function. It’s called mathemagical for a reason - you almost can’t believe that it works.
But mulDiv
works only with uint256
numbers. To go around this, I compute the signs and the absolute values.
Exp
Signature
function exp(int256 x) internal pure returns (int256 result);
Explanation
I use the following mathematical trick:
I already have an implementation for 2^x, and the second factor I can pre-compute offchain and use as a run-time constant. The binary logarithm of e is 1.442695040888963407.
Caveats:
- x must be less than 88.722839111672999628.
- For any x less than -41.446531673892822322, the result is zero.
Exp2
Signature
function exp2(int256 x) internal pure returns (int256 result);
Explanation
The key insight here is to write the binary exponential as a binary fraction expansion.
Any fractional number can be written as the sum of the integer part and the fractional part:
Then, using x as an exponent:
Writing f as a binary fraction:
Where f_i can be either 0 or 1. Then:
Note that:
I pre-compute these magic factors \sqrt[2^i]{2} off-chain use them as run-time constants. The i stands for the position of the bit in the fractional binary representation. E.g. if I were to match the very first bit, I would multiply by the square root of 2, i.e. ~1.4142.
Caveats:
- x must be lower than 128, because 2^128 doesn’t fit within the 128.128-bit fixed-point representation used internally in this function.
- When x is negative, I return the inverse of the binary exponential of the absolute value.
- For any x less than -59.794705707972522261, the result is zero.
Floor
Signature
function floor(int256 x) internal pure returns (int256 result);
Explanation
This function can be thought of an analogue of ceil
.
The first thing I do is check that x is greater than or equal to the value dubbed “min whole sd59x18”. This is min sd59x18 with all 18 decimals set to zero and the last digit incremented. This requirement is in place because the floor couldn’t otherwise be represented in Solidity.
Then, I calculate the remainder of x divided by the scale number 1e18, and:
- If the remainder is zero, return x.
- If the remainder is non-zero, I set the result to the difference between x and the remainder calculated above. And finally, if x is negative, I also subtract the scale number from the result.
Frac
Signature
function frac(int256 x) internal pure returns (int256 result);
Explanation
This is as simple as taking the modulo of x with respect to the scale number 1e18. The implementation is based on the odd function definition of the fractional part.
Gm
Signature
function gm(int256 x, int256 y) internal pure returns (int256 result);
Explanation
I first ensure that the product of the two inputs:
- Does not overflow.
- Is not negative.
Then, I calculate the square root of the product.
Inv
Signature
function inv(int256 x) internal pure returns (int256 result);
Explanation
I divide 1e36 (the square of the scale number) by the input.
Ln
Signature
function ln(int256 x) internal pure returns (int256 result);
Explanation
I perform some boundary checks and then I use this mathematical trick:
The first factor is the \log_2 function for which I already have an implementation, and second factor I pre-compute off-chain and I use as a run-time constant. The binary logarithm of e is 1.442695040888963407.
Caveats
- Doesn’t return exactly 1 for 2.718281828459045235, due to the lossy precision of the iterative approximation.
Log10
Signature
function log10(int256 x) internal pure returns (int256 result);
Explanation
Due to the lossy precision of the iterative approximation (explained below under the log2
section), the last decimal of the result, that is, the 18th decimal, may not be the same as when calculating the logarith with an advanced calculator like Wolfram.
This is an incovenience that can be glossed over in most cases, except when passing a power of ten as the input x
. In this case, the function should return a whole number. The algorithm goes something like this:
- If the number is a power of ten, return the power.
- If the number is not a power of ten, use this mathematical identity:
Checking that the input is a power of ten was easier said than done (I documented) my efforts toward finding the fastest possible implementation). The gist is that there are quite a few ways to verify that a number is power of ten, but most are achingly slow in Solidity. The most efficient one is to use assembly and the switch control operator and compare x against all hardcoded powers of ten that fit within the given type.
Log2
Signature
function log2(int256 x) internal pure returns (int256 result);
Explanation
I begin by looking at the sign of x:
- If it’s positive, do nothing.
- If it’s negative, take the inverse of x.
I can do this because of this mathematical identity:
Now I begin the iterative approximation algorithm. I will give a succinct description below, but do check the Wikipedia article I linked for a thorough explanation.
Let n be the integer part of the logarithm. Then, the fractional part of the logarithm, f, can be calculated like this:
Now, given the equality:
It follows that:
Thus, calculating the fractional part of a binary logarithm can be reduced to calculating the binary logarithm of a number between 1 (inclusive) and 2 (exclusive). To do this I use the following two rules:
The idea is to continuously apply the first rule until \frac{x}{2^n} crosses the value of two. Then, I apply the second rule and halve x.
Caveats
- The results are not perfectly accurate to the last decimal, due to the lossy precision of the iterative approximation.
Mul
Signature
function mul(int256 x, int256 y) internal pure returns (int256 result);
Explanation
Multiplying two fixed-point numbers is done by multiplying the two fixed-point numbers together, and then dividing the product by the scale number 1e18.
Just like in div
, I use Remco Bloemen’s mulDiv to fend off phantom overflow, and I compute the signs and the absolute values separately. But there’s two subtle things that I make differently here.
- I am not using
mulDiv
, but a variant that I wrote,mulDivFixedPoint
. The latter is twice more efficient than the former, because I replaced many runtime operations with constants pre-computed off-chain. The cool kids would call this constant folding. - I add 1 to the result if this property holds:
The goal with this is to increase the accuracy of the result in certain cases where Solidity would incorrectly truncate the result down. For example, without my optimization, a result of 6.6e-19 would be truncated to 0 instead of rounded up to 1 (recall that PRBMath works with maximum 18 decimals). For more details, see “listing 6” and text above it here.
If there’s anything that might be called innovative about PRBMath, it’s this mul
function. I basically took Remco’s function and refined it for fixed-point number multiplication.
Pow
Signature
function pow(int256 x, int256 y) internal pure returns (int256 result);
Explanation
This is an emergent function - it is an integration of other functions already defined in the library:
That is, exp2
, log2
and mul
.
Powu
Signature
function powu(int256 x, uint256 y) internal pure returns (int256 result);
Explanation
This is an alternative to pow
that works strictly with basic integer exponents, and it is slightly more precise.
I implemented this function using the famous algorithm exponentiation by squaring. Given the popularity of the algorithm, I’m sure there is already a much better explanation somewhere else on the Internet; I’ll nonetheless give a succinct description here.
Suppose n is the exponent. When n is odd, write x^n as:
And when n is even, write x^n as:
I begin by setting n = y and proceed with halfing n and applying one of the functions above, depending upon the parity of n.
Caveats
- This is the only function in PRBMath that doesn’t work only with fixed-point numbers. The second argument y is a basic unsigned integer.
- I assume that 0^0 is equal to 1.
Sqrt
Signature
function sqrt(int256 x) internal pure returns (int256 result);
Explanation
The idea here is to employ the Babylonian method and set the initial guess to the closest power of two that is higher than x.
I do seven iterations and I stop.
Gas Efficiency
You can find gas efficiency tables on GitHub. I won’t replicate them here to avoid having to make revisions in multiple locations.
Gas-wise, PRBMath is comparable to ABDKMath, which is reputed to be the fastest advanced math library. But given that PRBMath is faster than ABDKMath among 7 out of all 12 common functions, I think it’s safe to assert that PRBMath has replaced ABDKMath as the most efficient advanced math library for Solidity.
Testing and Security
There are 1,422 test cases and test coverage is around 96%, although in reality I think that the number is lower because solidity-coverage doesn’t cover unchecked
blocks. Nonetheless, the percentage should be pretty high.
To increase security guarantees, I plan to:
- Write fuzzing tests with Echidna.
- Get a security audit.
User Documentation
I don’t have a documentation website, but I went above and beyond to ensure that every single block of code is properly commented. In fact, I think that there are more commented lines than actual lines of code.
All functions are preceded by comments that adhere to the NatSpec Format, which makes the code clear, intuitive and sensible. There are even sections for requirements and caveats, some of which are replicated above.
You don’t need an undergrad in mathematics or computer science to understand the inner workings of PRBMath.
End Notes
PRBMath is released under Unlicense, but I don’t give any warranties and I won’t be liable for any loss, direct or indirect, through continued use.
Issues and PRs are welcome on GitHub. Questions I’d like to address the SCR community:
- Generally speaking, what do you think of this? Would you use PRBMath in your projects?
- Can you find any critical flaw anywhere?
- Do you know how to make any of my function implementations simpler, more efficient, or more elegant?
- What functions are missing? Why?