Deep Diving into PRBMath, a Library for Advanced Fixed-Point Math

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:

  1. PRBMathSD59x18
  2. PRBMathSD59x18Typed
  3. PRBMathUD60x18
  4. 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:

  1. If the remainder is zero, return x.
  2. 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:

e^x = 2^{x * \log_2{e}}

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:

x = n + f, 0 <= f < 1

Then, using x as an exponent:

y = 2^x = 2^{n + f} = 2^n * 2^f

Writing f as a binary fraction:

f = f1 * 2^{-1} + f2 * 2^{-2} + ...

Where f_i can be either 0 or 1. Then:

y = 2^n * 2^{f1*2^{-1}} * 2^{f2*2^{-2}} ...

Note that:

2^{f^i*2^{-1}} = (2^{2^{-i}})^{f_i} = \sqrt[2^i]{2^{f_i}} = \begin{cases} 1, & \text{if $f_i$ is 0} \\ \sqrt[2^i]{2}, & \text{if $f_i$ is 1} \end{cases}

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:

  1. If the remainder is zero, return x.
  2. 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:

  1. Does not overflow.
  2. 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:

\ln_x = \frac{\log_2{x}}{\log_2{e}}

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:

  1. If the number is a power of ten, return the power.
  2. If the number is not a power of ten, use this mathematical identity:
\log_{10}{x} = \frac{\log_2{x}}{\log_2{10}}

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:

  1. If it’s positive, do nothing.
  2. If it’s negative, take the inverse of x.

I can do this because of this mathematical identity:

\log_2{x} = -\log_2{\frac{1}{x}}

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:

f = \log_2{x} - n
f = \log_2{x} - log_2{2^n}
f = \frac{\log_2{x}}{2^n}

Now, given the equality:

2^n <= x < 2^{n+1}

It follows that:

1 <= \frac{x}{2^n} < 2

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:

\log_2{x} = \frac{\log_2{x^2}}{2}
\log_2{x} = 1 + \log_2{\frac{x}{2}}

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.

  1. 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.
  2. I add 1 to the result if this property holds:
(x * y) \bmod scale >= \frac{scale}{2}

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:

x^y = 2^{log2(x)*y}

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:

x*(x^2)^{\frac{n-1}{2}}

And when n is even, write x^n as:

(x^2)^{\frac{n}{2}}

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.

x_{n+1} = \frac{x_n+\frac{S}{x_n}}{2}

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:

  1. Write fuzzing tests with Echidna.
  2. 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:

  1. Generally speaking, what do you think of this? Would you use PRBMath in your projects?
  2. Can you find any critical flaw anywhere?
  3. Do you know how to make any of my function implementations simpler, more efficient, or more elegant?
  4. What functions are missing? Why?
15 Likes

Thanks for this great library and your explanations behind each function. Given all the benefits, if someone was building a financial application in solidity what reasons could they have for not using such a library?

5 Likes

If someone was building a financial application in solidity what reasons could they have for not using such a library?

If you’re using Solidity v0.7 or below, you can’t use PRBMath. But for v0.8, indeed, this is a solid choice for a fixed-point math library in a Solidity project.

There is one edge case though. If you care about every single gas unit, and you happen to only need to use either powu or sqrt, it might be cheaper to use ABDKMath. Albeit you will have to work with 128.128-bit numbers, it has faster implementations for these specific functions.

6 Likes

@paulrberg Thanks for this post. This is a joy to read.

A couple of points to clarify:

  1. How close is the current implementation from production usage? Are there any projects using it?
  2. Looking at coveralls report on Github, while the overall test coverage is high (line coverage at is 95% and branch coverage at almost 97%), I noticed that for contracts/PRBMath.sol, the overall coverage is at 83%. Any particular reason for that drop in comparison to the other files?
  3. As you pointed out, solidity coverage is not counting unchecked blocks, as their content does not appear green in the coveralls report. Is your test suite indeed stressing those paths?
  4. Was the entire code written from scratch (seems like), or was it partially cloned from another library? If so, having this sort of traceability documented in the code would help towards a future audit.

Please let me know of your thoughts.

6 Likes

GitHub reports that there are currently 21 projects using PRBMath:

GitHub Users

What portion of these are serious projects, and what are short-lived hackathon experiments, I can’t know for sure. What I do know is that we at Hifi will use the library in our v1 release. Also the Pods team, who I met in person, told me that they are looking into PRBMath for their next release as well.

Note that the core PRBMath.sol file is not meant to be used directly by end-users. I wrote it with the specific purpose of avoiding duplicate code.

Coverage is ~83% because the mulDivSigned function is not tested. But this function is not used in SD59x18, SD59x18Typed, UD60x18 or UD60x18Typed. I added it an “alpha” feature that can be used by experienced users who know what they are doing (none of the instructions in the README contain an import for PRBMath.sol, so in theory users should not know about it).

It certainly does. The vast majority of my usage of unchecked is “exhaustive”, meaning that the entire function is wrapped within an unchecked block.

Side note: this is when I realised that it would be great if Solidity added an “unchecked” keyword that we could add as a modifier to functions that should have disabled the usual overflow/ underflow checks.

From scratch, but I need to expound on this.

While some of the features are entirely novel, like the avg function or the different mulDiv variants (mulDivFixedPoint, mulDivSigned), I sourced the mathematical logic behind some functions (e.g. exp2, pow) from Mikhail Vladimirov’s Math in Solidity series. Mikhail is the author of ABDKMath, but his article series does not share snippets from the library; it just explains the math.

This is conjecture on my part, but I suspect that he also sourced a part of his implementations from Wikipedia or other places. Some mathematical functions can most efficiently be implemented with certain known algorithms, which are timeless. E.g. the Babylonian method for the square root.

5 Likes

@paulrberg Thanks for the clarification and effort put in the explanation. All is clear :)

Still, if you allow me, I would like to follow-up on the first question you had in your post (questions 2 & 3 are sort of by-products):

I think this is a great project and the value is there. Not only it brings more functionality than ABDK, but it is also more gas-efficient for certain operations. However, ABDK’s library has been extensively used by now, and since ABDK is itself an auditing company, there is some level of assumed trust. Hence, as you acknowledge yourself, the next step is to make sure your code is secure as best as you can. Consider setting up a crowd-funding arrangement to get an audit for this library. Without such an audit in place, adoption could be hindered.

Also, you mentioned using Echidna to fuzz test and bring more security guarantees/confidence that the implementation is not buggy:

Note: Scribble could potentially be an alternative; @maurelian wrote a nice post on SCRF about it. In a nutshell, “Scribble is a specification language and runtime verification tool that translates high-level specifications into Solidity code” (ref: Scribble | Consensys Diligence). Following their website, “[…] after writing properties, developers can use tools such as Diligence Fuzzing to automatically test smart contracts and ensure all is working as planned!”. Their fuzzing as a service tool is currently being released as an early access style, but should be something to be considered.

Hope that helps. Once again thanks for your post and welcome to SCRF :slight_smile:

2 Likes

I thank you! This was fantastic feedback.

Yes, I actually have an idea in mind for how to run the crowdfund. I agree that an audit would go a long way.

Interesting, I didn’t know about it, thanks. I will check it out!

4 Likes

I just saw this Solidity rounding error in Mean Finance. @paulrberg do you avoid this kind of bugs?
More generally, how do you guarantee that a rounding error somewhere won’t snowball into a bigger error somewhere else?

4 Likes

I did all I could to make the library as precise as possible.

Firstly, take a look at mulDivFixedPoint, which contains some hardcore assembly code. What that does is to implement the logic explained in listing 6 of this journal article, i.e. to increase the accuracy of the mulDiv result in cases where Solidity would otherwise truncate the result down. Without my adjustment, a value like 6.6e-19 would be rounded down 0, instead of being rounded up to 1, which is numerically closer.

Secondly, most PRBMath functions have a “Caveats” section in their NaSpec that explain what shortcomings they have. Here’s an example for the ln (natural logarithm) function:

This doesn’t return exactly 1 for 2718281828459045235, for that we would need more fine-grained precision.

Now, about this:

There’s so much precision you can get in Solidity.

Any Solidity library developer who guarantees that there aren’t any rounding errors is either a fool or a fraud. The onus is on the library user to ensure that their protocol is safe to use.

6 Likes

Thank you for the library and the writeup, Paul!

I am to create a vault that gives out rewards based off of the geometric sequence. To do this, i need to perform operations like 0.9993^10000
I am having trouble making this happen with your library for base values less then 1. Here is an example that tries to calculate 0.9993^3, which reverts. I also tried with powu, which returns 0.

  function doTheMath() external pure returns (uint256 result) {
    uint256 x = 9993*10**14;
    uint256 y = 3*10**18;
    result = (x).pow(y);
  }
2 Likes

Hi @djrthree, I’m glad that you like it!

The best channel to ask for support is the Discussions section on GitHub.

I’d like to not offer any answer to your question here, because the behavior of the pow function might change in the future for sub-unitary bases.

3 Likes

Update: I have recently released PRBMath v3, which is a significant upgrade with many breaking changes.

What’s new:

  • User defined value types SD59x18 and UD60x18
  • Free functions
  • Migration to Foundry
  • Various optimizations

You can read the full changelog for V3 in the release notes on GitHub, and see more code snippet examples in the README.

The mathematical explanations given above should still stand, but you should note that now all PRBMath functions use the value types SD59x18 and UD60x18 instead of int256 and uint256, respectively.

1 Like