Gibt es eine effiziente Möglichkeit, die Potenzierung eines Bruchs und einer ganzen Zahl zu berechnen?

Beachten Sie die folgende JavaScript-Operation:

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

Dabei rist ein Bruch (z. B. 1.01) und neine ganze Zahl. Gibt es eine effiziente Möglichkeit, diese Operation effizient auf Ethereum zu emulieren? Das heißt, eine Funktion:

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

Welche gibt ähnliche Werte wie die erstere zurück? Eine solche Funktion hätte viele Anwendungsfälle; beispielsweise die Berechnung von Zinseszinssätzen.

Sie benötigen wahrscheinlich die ABDK MAth 64.64- Bibliothek. Es arbeitet mit binären Festkommazahlen (64 Binärziffern nach dem Punkt) und verfügt über alle grundlegenden mathematischen Operationen einschließlich Potenzoperationen.

Antworten (3)

Näherungslösung unter Verwendung der Binomialentwicklung

Das Folgende ist eine anständige, kostengünstige Annäherung:

// 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;
}

Diese Funktion berechnet k * (1+1/q) ^ N. Um also beispielsweise zu berechnen 2500 * 1.01 ^ 137, könnten wir verwenden fracExp(2500, 100, 137, 8). Dies gibt 9769aus, was dem realen Wert sehr nahe kommt 9771.657061221898.


Erläuterung

1. Darstellung rationaler Zahlen mit Fixpunkten

Das Hauptproblem ist, dass Ethereum keine Gleitkommazahlen hat. Der übliche Ansatz zur Darstellung reeller Zahlen ist die Verwendung von Fixpunkten. Das heißt, wir multiplizieren alles mit einer großen Zahl und nehmen an, dass es einen impliziten Punkt gibt. Beispielsweise 1.002wird dargestellt als 1002. Für die meisten arithmetischen Operationen brauchen wir eigentlich keine Gleitkommazahlen; Beispiele:

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

2. Festkommapotenzierung

Für die Potenzierung funktioniert dieser Ansatz jedoch nicht so gut. Um es zu implementieren, müssten wir den Zähler und den Nenner potenzieren und dann dividieren. Beispiel: 1.1^75 == (11/10)^75 == (11^75)/(10^75) == 1271. Das Problem ist, dass wir, um zum Ergebnis zu gelangen, den Zwischenwert berechnen 11^75mussten, der größer als 2^256ist, was zu einem Integer-Überlauf führte . Beide Operanden und die Ergebnisse könnten in ein 16-Bit-Uint passen, aber wir überlaufen die 256-Bit-Wortgröße!

Eine Möglichkeit, den Überlauf zu vermeiden, besteht darin, die Potenzierung als Schleife zu implementieren. Bei jedem Durchlauf multiplizieren Sie die Zahl und normalisieren sie neu, damit sie nicht zu groß wird. Das Problem bei diesem Ansatz ist, dass er exponentiell linear ist und daher unerschwinglich teuer sein kann. Es gibt effizientere exp-Implementierungen wie Potenzierung durch Quadrieren ; Ich bin mir jedoch nicht sicher, ob dieser Ansatz die Normalisierung der Zwischenbedingungen ermöglicht.

3. Vermeidung großer Zwischenwerte durch Binomialentwicklung

Schauen wir uns noch einmal ein Beispiel für die Art von Berechnung an, die wir durchzuführen versuchen:

1.01 ^ 100

Wir können dies umschreiben als:

(0.01 + 1) ^ 100

Wendet man die Binomialentwicklung an, wird das zu:

  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

Da 1^Nimmer 1 ist, ist dies dasselbe wie:

  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

Beachten Sie, dass jede neue Zeile zunehmend kleinere Begriffe enthält. Aus diesem Grund können wir einen guten Teil dieser Linien eliminieren, ohne das Ergebnis erheblich zu beeinträchtigen. Diese Lösung gibt uns eine vernünftige Annäherung. Die hier gepostete Lösung berechnet die persten Zeilen dieser Reihe.

"Wenn p == n, ist die Genauigkeit absolut" ist dies NICHT wahr. Der Code approximiert das Binomial über eine unendliche Reihenentwicklung. Je höher p ist, desto genauer ist die Schätzung, aber keine endliche Zahl für p ergibt "absolute" Genauigkeit.

Aus diesem Vertrag: 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;
    }
  }
}

Dies berechnet x^y, wobei x eine 128,128-Bit-Festkommazahl in Wut [0..1) und y eine 256-Bit-Ganzzahl ohne Vorzeichen ist.

Es sollte nicht allzu schwer sein, diesen Code zu ändern, um x in einem größeren Bereich zuzulassen.

Lösung

Meine erweiterte Mathematikbibliothek PRBMath bietet das, wonach Sie suchen, über die powFunktion, die die Potenzierung durch Quadrierungsalgorithmus implementiert . Ich werde meine Implementierung hier für die Nachwelt einfügen, aber sehen Sie sich das verlinkte Repo für den aktuellsten Code an.

/// @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

Die obige Funktion geht von der vorzeichenbehafteten 59,18-Dezimal-Festkommadarstellung (kurz gesagt als SD59x18 geschrieben) aus, bei der es sich um Ganzzahlen handelt, die Festkomma emulieren, indem sie die letzten 18 Ziffern als die ersten 18 Dezimalstellen betrachten. Zum Beispiel ewird geschrieben als 2718281828459045235.