-
Notifications
You must be signed in to change notification settings - Fork 48
/
pow_approx.go
98 lines (85 loc) · 2.87 KB
/
pow_approx.go
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
package types
import (
"fmt"
sdk "github.com/cosmos/cosmos-sdk/types"
)
// Contract: 0 < base <= 2
// 0 <= exp < 1.
func PowApprox(base sdk.Dec, exp sdk.Dec, precision sdk.Dec) sdk.Dec {
if !base.IsPositive() {
panic(fmt.Errorf("base must be greater than 0"))
}
if exp.IsZero() {
return sdk.OneDec()
}
// Common case optimization
// Optimize for it being equal to one-half
if exp.Equal(one_half) {
output, err := base.ApproxSqrt()
if err != nil {
panic(err)
}
return output
}
// TODO: Make an approx-equal function, and then check if exp * 3 = 1, and do a check accordingly
// We compute this via taking the maclaurin series of (1 + x)^a
// where x = base - 1.
// The maclaurin series of (1 + x)^a = sum_{k=0}^{infty} binom(a, k) x^k
// Binom(a, k) takes the natural continuation on the first parameter, namely that
// Binom(a, k) = N/D, where D = k!, and N = a(a-1)(a-2)...(a-k+1)
// Next we show that the absolute value of each term is less than the last term.
// Note that the change in term n's value vs term n + 1 is a multiplicative factor of
// v_n = x(a - n) / (n+1)
// So if |v_n| < 1, we know that each term has a lesser impact on the result than the last.
// For our bounds on |x| < 1, |a| < 1,
// it suffices to see for what n is |v_n| < 1,
// in the worst parameterization of x = 1, a = -1.
// v_n = |(-1 + epsilon - n) / (n+1)|
// So |v_n| is always less than 1, as n ranges over the integers.
//
// Note that term_n of the expansion is 1 * prod_{i=0}^{n-1} v_i
// The error if we stop the expansion at term_n is:
// error_n = sum_{k=n+1}^{infty} term_k
// At this point we further restrict a >= 0, so 0 <= a < 1.
// Now we take the _INCORRECT_ assumption that if term_n < p, then
// error_n < p.
// This assumption is obviously wrong.
// However our usages of this function don't use the full domain.
// With a > 0, |x| << 1, and p sufficiently low, perhaps this actually is true.
// TODO: Check with our parameterization
// TODO: If theres a bug, balancer is also wrong here :thonk:
base = base.Clone()
x, xneg := AbsDifferenceWithSign(base, one)
term := sdk.OneDec()
sum := sdk.OneDec()
negative := false
a := exp.Clone()
bigK := sdk.NewDec(0)
// TODO: Document this computation via taylor expansion
for i := int64(1); term.GTE(precision); i++ {
// At each iteration, we need two values, i and i-1.
// To avoid expensive big.Int allocation, we reuse bigK variable.
// On this line, bigK == i-1.
c, cneg := AbsDifferenceWithSign(a, bigK)
// On this line, bigK == i.
bigK.Set(sdk.NewDec(i)) // TODO: O(n) bigint allocation happens
term.MulMut(c).MulMut(x).QuoMut(bigK)
// a is mutated on absDifferenceWithSign, reset
a.Set(exp)
if term.IsZero() {
break
}
if xneg {
negative = !negative
}
if cneg {
negative = !negative
}
if negative {
sum.SubMut(term)
} else {
sum.AddMut(term)
}
}
return sum
}