是否有一种有效的方法来近似 (a / b)^n,其中 a、b 和 n 是无符号整数?
Is there an efficient way to approximate (a / b)^n where a, b, and n are unsigned integers?
Exponentiation by squaring是一种快速计算an的算法,其中a
和n
是有符号整数。 (它在 O(log n) 乘法中这样做)。
是否有类似的算法,而是计算 (a / b)n,其中 a
、b
和 n
都是无符号整数吗?显而易见的方法(即计算 an / bn)的问题在于它会 return 由于整数导致的错误结果中间值溢出。
宿主语言中没有浮点数,只有整数。
我可以提供一个大概的答案。
如果您希望 (a/b)^n 的值具有出色的准确性,其中 a、b 和 n 是无符号整数并且您没有可用的浮点运算--请使用扩展精度整数计算求出a^n和b^n,然后将两者相除。
某些语言(例如 Python)内置了扩展精度整数算法。如果您的语言没有,请寻找实现它的程序包。如果你做不到,那就自己做一个包。这并不难——这样的软件包是我当时第二学期计算机科学 class 的作业。乘法和幂相当简单;最困难的部分是除法,即使你只想要商和余数。但是 "most difficult" 并不意味着 "very difficult",您可能可以做到。第二个必须困难的例程是将扩展整数打印为十进制格式。
基本思想是将每个整数存储在数组或常规无符号整数列表中,其中整数是 "digit" 具有大基数的算术。您希望能够处理任意两位数字的乘积,因此如果您的机器有 32 位整数而您无法处理 64 位整数,请存储每个 16 位的 "digits"。 "digit" 越大,计算速度越快。如果您的计算很少,并且经常打印到小数,请使用 10 的幂,例如每个 "digit".
使用 10000
询问您是否需要更多详细信息。
以防万一有人正在寻找常数-space 解决方案,我已经用二项式展开解决了这个问题,这是一个不错的近似值。我正在使用以下代码:
// 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.
上发布了更周到的解释
这是一个基于费曼对数算法的定点 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
}
Exponentiation by squaring是一种快速计算an的算法,其中a
和n
是有符号整数。 (它在 O(log n) 乘法中这样做)。
是否有类似的算法,而是计算 (a / b)n,其中 a
、b
和 n
都是无符号整数吗?显而易见的方法(即计算 an / bn)的问题在于它会 return 由于整数导致的错误结果中间值溢出。
宿主语言中没有浮点数,只有整数。
我可以提供一个大概的答案。
如果您希望 (a/b)^n 的值具有出色的准确性,其中 a、b 和 n 是无符号整数并且您没有可用的浮点运算--请使用扩展精度整数计算求出a^n和b^n,然后将两者相除。
某些语言(例如 Python)内置了扩展精度整数算法。如果您的语言没有,请寻找实现它的程序包。如果你做不到,那就自己做一个包。这并不难——这样的软件包是我当时第二学期计算机科学 class 的作业。乘法和幂相当简单;最困难的部分是除法,即使你只想要商和余数。但是 "most difficult" 并不意味着 "very difficult",您可能可以做到。第二个必须困难的例程是将扩展整数打印为十进制格式。
基本思想是将每个整数存储在数组或常规无符号整数列表中,其中整数是 "digit" 具有大基数的算术。您希望能够处理任意两位数字的乘积,因此如果您的机器有 32 位整数而您无法处理 64 位整数,请存储每个 16 位的 "digits"。 "digit" 越大,计算速度越快。如果您的计算很少,并且经常打印到小数,请使用 10 的幂,例如每个 "digit".
使用 10000询问您是否需要更多详细信息。
以防万一有人正在寻找常数-space 解决方案,我已经用二项式展开解决了这个问题,这是一个不错的近似值。我正在使用以下代码:
// 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.
这是一个基于费曼对数算法的定点 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
}