有没有一种有效的方法来近似 (a / b)^n，其中 a、b 和 n 是无符号整数？

Question

Exponentiation by squaring 是一种快速计算aⁿ 的算法，其中a 和n 是有符号整数。 （它在 O(log n) 乘法中这样做）。

是否有类似的算法，而是计算 (a / b)ⁿ，其中 a、b 和 n 都是无符号整数？显而易见的方法（即计算 aⁿ / bⁿ）的问题在于，由于中间值的整数溢出，它会返回错误的结果。

我在宿主语言中没有浮点数，只有整数。

我可以给出一个大概的答案。

Answer 1

如果您希望 (a/b)^n 的值具有出色的精度，其中 a、b 和 n 是无符号整数，并且您没有可用的浮点运算 - 使用扩展精度整数计算来找到^n 和 b^n，然后将两者相除。

某些语言（例如 Python）内置了扩展精度整数运算。如果您的语言没有它，请寻找实现它的包。如果你做不到，那就自己做一个包。这并不难——这样一个包是我当时第二学期计算机科学课上的一个作业。乘法和幂相当简单。最困难的部分是除法，即使你只想要商和余数。但是“最困难”并不意味着“非常困难”，您可能可以做到。第二个必须困难的例程是将扩展整数打印为十进制格式。

基本思想是将每个整数存储在一个数组或常规无符号整数列表中，其中整数是具有大基数的算术“数字”。您希望能够处理任何两位数的乘积，因此如果您的机器有 32 位整数并且您无法处理 64 位整数，请存储每个 16 位的“数字”。 “数字”越大，计算速度越快。如果您的计算量很少并且经常打印到小数点，请为每个“数字”使用 10 的幂，例如 10000。

询问您是否需要更多详细信息。

Answer 2

这是一个基于 Feynman 对数算法的定点 pow 实现。它很快而且有点脏； C 库倾向于使用多项式近似，但这种方法更复杂，我不确定它转换为定点的效果如何。

// powFraction approximates (a/b)**n.
func powFraction(a uint64, b uint64, n uint64) uint64 {
    if a == 0 || b == 0 || a < b {
        panic("powFraction")
    }
    return expFixed((logFixed(a) - logFixed(b)) * n)
}

// logFixed approximates 2**58 * log2(x). [Feynman]
func logFixed(x uint64) uint64 {
    if x == 0 {
        panic("logFixed")
    }
    // Normalize x into [2**63, 2**64).
    n := numberOfLeadingZeros(x)
    x <<= n
    p := uint64(1 << 63)
    y := uint64(0)
    for k := uint(1); k <= 63; k++ {
        // Warning: if q > x-p, then p + q may overflow.
        if q := p >> k; q <= x-p {
            p += q
            y += table[k-1]
        }
    }
    return uint64(63-n)<<58 + y>>6
}

// expFixed approximately inverts logFixed.
func expFixed(y uint64) uint64 {
    n := 63 - uint(y>>58)
    y <<= 6
    p := uint64(1 << 63)
    for k := uint(1); k <= 63; k++ {
        if z := table[k-1]; z <= y {
            p += p >> k
            y -= z
        }
    }
    return p >> n
}

// numberOfLeadingZeros returns the number of leading zeros in the word x.
// [Hacker's Delight]
func numberOfLeadingZeros(x uint64) uint {
    n := uint(64)
    if y := x >> 32; y != 0 {
        x = y
        n = 32
    }
    if y := x >> 16; y != 0 {
        x = y
        n -= 16
    }
    if y := x >> 8; y != 0 {
        x = y
        n -= 8
    }
    if y := x >> 4; y != 0 {
        x = y
        n -= 4
    }
    if y := x >> 2; y != 0 {
        x = y
        n -= 2
    }
    if x>>1 != 0 {
        return n - 2
    }
    return n - uint(x)
}

// table[k-1] approximates 2**64 * log2(1 + 2**-k). [MPFR]
var table = [...]uint64{
    10790653543520307104, // 1
    5938525176524057593,  // 2
    3134563013331062591,  // 3
    1613404648504497789,  // 4
    818926958183105433,   // 5
    412613322424486499,   // 6
    207106307442936368,   // 7
    103754619509458805,   // 8
    51927872466823974,    // 9
    25976601570169168,    // 10
    12991470209511302,    // 11
    6496527847636937,     // 12
    3248462157916594,     // 13
    1624280643531991,     // 14
    812152713665686,      // 15
    406079454902306,      // 16
    203040501980337,      // 17
    101520444623942,      // 18
    50760270720599,       // 19
    25380147462480,       // 20
    12690076756788,       // 21
    6345039134781,        // 22
    3172519756487,        // 23
    1586259925518,        // 24
    793129974578,         // 25
    396564990243,         // 26
    198282495860,         // 27
    99141248115,          // 28
    49570624104,          // 29
    24785312063,          // 30
    12392656035,          // 31
    6196328018,           // 32
    3098164009,           // 33
    1549082005,           // 34
    774541002,            // 35
    387270501,            // 36
    193635251,            // 37
    96817625,             // 38
    48408813,             // 39
    24204406,             // 40
    12102203,             // 41
    6051102,              // 42
    3025551,              // 43
    1512775,              // 44
    756388,               // 45
    378194,               // 46
    189097,               // 47
    94548,                // 48
    47274,                // 49
    23637,                // 50
    11819,                // 51
    5909,                 // 52
    2955,                 // 53
    1477,                 // 54
    739,                  // 55
    369,                  // 56
    185,                  // 57
    92,                   // 58
    46,                   // 59
    23,                   // 60
    12,                   // 61
    6,                    // 62
    3,                    // 63
}

Answer 3

以防万一有人在寻找恒定空间解决方案，我已经通过二项式展开解决了这个问题，这是一个不错的近似值。我正在使用以下代码：

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

这只是计算(1 + r)^N 的二项式展开的p 第一项，其中r 是一个小的正实数。我在Ethereum Stack Exchange发布了一个更周到的解释。