13

Mind the following JavaScript operation:

function exp(n){
    return Math.floor(Math.pow(1.01, n));
}

Where r is a fraction (such as 1.01) and n is an integer. Is there any efficient way to emulate that operation efficiently on Ethereum? That is, a function:

function exp(uint n) returns (uint) {
    return ...
}

Which returns similar values to the former? Such a function would have many use cases; for example, computing compound interest rates.

1
  • You probably need ABDK MAth 64.64 library. It operates with binary fixed point numbers (64 binary digits after dot) and has all basic math operations including power operation. – Mikhail Vladimirov Apr 11 '19 at 11:05
16

Approximate solution using the binomial expansion

The following is a decent, low-cost approximation:

// Computes `k * (1+1/q) ^ N`, with precision `p`. The higher
// the precision, the higher the gas cost. It should be
// something around the log of `n`. When `p == n`, the
// precision is absolute (sans possible integer overflows). <edit: NOT true, see comments>
// Much smaller values are sufficient to get a great approximation.
function fracExp(uint k, uint q, uint n, uint p) returns (uint) {
  uint s = 0;
  uint N = 1;
  uint B = 1;
  for (uint i = 0; i < p; ++i){
    s += k * N / B / (q**i);
    N  = N * (n-i);
    B  = B * (i+1);
  }
  return s;
}

This function computes k * (1+1/q) ^ N. So, for example, to compute 2500 * 1.01 ^ 137, we could use fracExp(2500, 100, 137, 8). This outputs 9769, which is very close to the real value, 9771.657061221898.


Explanation

1. Representing rational numbers with fixed points

The main issue is Ethereum doesn't have floating points. The usual approach to represent real numbers is by using fixed points. That means we multiply everything by a large number and assume there is an implicit dot. For example, 1.002 is represented as 1002. For most arithmetic operations we don't really need floats; examples:

  • Addition: 1.1 + 1.2 == 11/10 + 12/10 == 23/10 == 2.3
  • Multiplication: 1.1 * 2 == 11/10 * 2 = 22/10 == 2.2

2. Fixed point exponentiation

For exponentiation, though, that approach doesn't work so well. To implement it, we'd need to exponentiate the numerator and the denominator, and then divide. Example: 1.1^75 == (11/10)^75 == (11^75)/(10^75) == 1271. The problem is that, to get to the result, we needed to compute the intermediate value 11^75, which is larger than 2^256, resulting in an integer overflow. Both operands and the results could fit in a 16-bit uint, yet we overflow the 256-bit word size!

A way to avoid the overflow is to implement exponentiation as a loop; at each pass, you multiply the number and re-normalize it so it doesn't get too big. The issue with that approach is that it is linear on the exponent, so it might be prohibitely expensive. There are more efficient exp implementations such as exponentiation by squaring; but I'm not sure that approach allows normalizing the intermediate terms.

3. Avoiding big intermediate values with binomial expansion

Let's look again at an example of the kind computation we're trying to perform:

1.01 ^ 100

We can rewrite this as:

(0.01 + 1) ^ 100

Applying the binomial expansion, that becomes:

  1                  * 1 ^ 100 * 0.1 ^ 0
+ 100                * 1 ^  99 * 0.1 ^ 1
+ 100*99/2           * 1 ^  98 * 0.1 ^ 2
+ 100*99*98/2/3      * 1 ^  97 * 0.1 ^ 3
+ 100*99*98*97/2/3/4 * 1 ^  96 * 0.1 ^ 4
...
+ 1                  * 1 ^   0 * 0.1 ^ 100

Since 1^N is always 1, this is the same as:

  1                  * 0.1 ^ 0
+ 100                * 0.1 ^ 1
+ 100*99/2           * 0.1 ^ 2
+ 100*99*98/2/3      * 0.1 ^ 3
+ 100*99*98*97/2/3/4 * 0.1 ^ 4
...
+ 1                  * 0.1 ^ 100

Notice that each new line contains progressively smaller terms. Due to that, we're able to eliminate a good chunk of those lines without affecting the result considerably. That solution gives us a reasonable approximation. The solution posted here computes the p first lines of that series.

1
  • 2
    "When p == n, the precision is absolute" this is NOT true. The code approximates the binomial via an infinite series expansion. The higher p is the more accurate the estimation, but no finite number for p will give "absolute" precision. – rohanp Jan 3 '18 at 22:45
1

From this contract: https://etherscan.io/address/0x21bceeef718a0928c2cc1f1d980bab5993f837ba#code

uint256 constant TWO_128 = 0x100000000000000000000000000000000; // 2^128
uint256 constant TWO_127 = 0x80000000000000000000000000000000; // 2^127

/**
 * Multiply _a by _b / 2^128.  Parameter _a should be less than or equal to
 * 2^128 and parameter _b should be less than 2^128.
 *
 * @param _a left argument
 * @param _b right argument
 * @return _a * _b / 2^128
 */
function mul (uint256 _a, uint256 _b)
internal constant returns (uint256 _result) {
  if (_a > TWO_128) throw;
  if (_b >= TWO_128) throw;
  return (_a * _b + TWO_127) >> 128;
}

/**
 * Calculate (_a / 2^128)^_b * 2^128.  Parameter _a should be less than 2^128.
 *
 * @param _a left argument
 * @param _b right argument
 * @return (_a / 2^128)^_b * 2^128
 */
function pow (uint256 _a, uint256 _b)
internal constant returns (uint256 _result) {
  if (_a >= TWO_128) throw;

  _result = TWO_128;
  while (_b > 0) {
    if (_b & 1 == 0) {
      _a = mul (_a, _a);
      _b >>= 1;
    } else {
      _result = mul (_result, _a);
      _b -= 1;
    }
  }
}

This calculates x^y where x is fixed point 128.128 bit number in rage [0..1) and y is 256-bit unsigned integer.

Should be not too hard to modify this code to allow x in wider range.

0
0

Solution

My advanced math library PRBMath offers what you're looking for via the pow function, which implements the exponentiation by squaring algorithm. I'm gonna paste my implementation here, for posterity, but check out the linked repo for the most up-to-date code.

/// @dev How many trailing decimals can be represented.
uint256 internal constant SCALE = 1e18;

/// @dev Largest power of two divisor of SCALE.
uint256 internal constant SCALE_LPOTD = 262144;

/// @dev SCALE inverted mod 2^256.
uint256 internal constant SCALE_INVERSE = 78156646155174841979727994598816262306175212592076161876661508869554232690281;

/// @notice Calculates floor(x*y÷1e18) with full precision.
///
/// @dev Variant of "mulDiv" with constant folding, i.e. in which the denominator is always 1e18. Before returning the
/// final result, we add 1 if (x * y) % SCALE >= HALF_SCALE. Without this, 6.6e-19 would be truncated to 0 instead of
/// being rounded to 1e-18.  See "Listing 6" and text above it at https://accu.org/index.php/journals/1717.
///
/// Requirements:
/// - The result must fit within uint256.
///
/// Caveats:
/// - The body is purposely left uncommented; see the NatSpec comments in "PRBMathCommon.mulDiv" to understand how this works.
/// - It is assumed that the result can never be type(uint256).max when x and y solve the following two queations:
///     1) x * y = type(uint256).max * SCALE
///     2) (x * y) % SCALE >= SCALE / 2
///
/// @param x The multiplicand as an unsigned 60.18-decimal fixed-point number.
/// @param y The multiplier as an unsigned 60.18-decimal fixed-point number.
/// @return result The result as an unsigned 60.18-decimal fixed-point number.
function mulDivFixedPoint(uint256 x, uint256 y) internal pure returns (uint256 result) {
    uint256 prod0;
    uint256 prod1;
    assembly {
        let mm := mulmod(x, y, not(0))
        prod0 := mul(x, y)
        prod1 := sub(sub(mm, prod0), lt(mm, prod0))
    }

    uint256 remainder;
    uint256 roundUpUnit;
    assembly {
        remainder := mulmod(x, y, SCALE)
        roundUpUnit := gt(remainder, 499999999999999999)
    }

    if (prod1 == 0) {
        unchecked {
            result = (prod0 / SCALE) + roundUpUnit;
            return result;
        }
    }

    require(SCALE > prod1);

    assembly {
        result := add(
            mul(
                or(
                    div(sub(prod0, remainder), SCALE_LPOTD),
                    mul(sub(prod1, gt(remainder, prod0)), add(div(sub(0, SCALE_LPOTD), SCALE_LPOTD), 1))
                ),
                SCALE_INVERSE
            ),
            roundUpUnit
        )
    }
}

/// @notice Raises x (unsigned 60.18-decimal fixed-point number) to the power of y (basic unsigned integer) using the
/// famous algorithm "exponentiation by squaring".
///
/// @dev See https://en.wikipedia.org/wiki/Exponentiation_by_squaring
///
/// Requirements:
/// - The result must fit within MAX_UD60x18.
///
/// Caveats:
/// - All from "mul".
/// - Assumes 0^0 is 1.
///
/// @param x The base as an unsigned 60.18-decimal fixed-point number.
/// @param y The exponent as an uint256.
/// @return result The result as an unsigned 60.18-decimal fixed-point number.
function pow(uint256 x, uint256 y) internal pure returns (uint256 result) {
    // Calculate the first iteration of the loop in advance.
    result = y & 1 > 0 ? x : SCALE;

    // Equivalent to "for(y /= 2; y > 0; y /= 2)" but faster.
    for (y >>= 1; y > 0; y >>= 1) {
        x = mulDivFixedPoint(x, x);

        // Equivalent to "y % 2 == 1" but faster.
        if (y & 1 > 0) {
            result = mulDivFixedPoint(result, x);
        }
    }
}

Nota Bene

The function above assumes the signed 59.18-decimal fixed-point representation (succinctly written as SD59x18), which are integer numbers that emulate fixed-point by considering the last 18 digits the first 18 decimals. For example, e is written as 2718281828459045235.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.