-
Notifications
You must be signed in to change notification settings - Fork 18
/
AnalyticMath.sol
178 lines (161 loc) · 10.7 KB
/
AnalyticMath.sol
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
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
// SPDX-License-Identifier: SEE LICENSE IN LICENSE
pragma solidity 0.8.28;
import "./IntegralMath.sol";
library AnalyticMath {
// Auto-generated via 'PrintAnalyticMathConstants.py'
uint256 internal constant FIXED_1 = 0x0080000000000000000000000000000000;
uint256 internal constant LN2_MIN = 0x0058b90bfbe8e7bcd5e4f1d9cc01f97b57;
uint256 internal constant LN2_MAX = 0x0058b90bfbe8e7bcd5e4f1d9cc01f97b58;
uint256 internal constant LOG_MID = 0x015bf0a8b1457695355fb8ac404e7a79e4;
uint256 internal constant EXP_MID = 0x0400000000000000000000000000000000;
uint256 internal constant EXP_MAX = 0x2cb53f09f05cc627c85ddebfccfeb72758;
/**
* @dev Compute (a / b) ^ (c / d)
*/
function pow(uint256 a, uint256 b, uint256 c, uint256 d) internal pure returns (uint256, uint256) { unchecked {
if (b == 0 || d == 0)
revert("division by zero");
if (a == 0 || c == 0)
return (a ** c, 1);
if (a > b)
return (fixedExp(IntegralMath.mulDivF(fixedLog(IntegralMath.mulDivF(FIXED_1, a, b)), c, d)), FIXED_1);
if (b > a)
return (FIXED_1, fixedExp(IntegralMath.mulDivF(fixedLog(IntegralMath.mulDivF(FIXED_1, b, a)), c, d)));
return (1, 1);
}}
/**
* @dev Compute log(a / b)
*/
function log(uint256 a, uint256 b) internal pure returns (uint256, uint256) { unchecked {
return (fixedLog(IntegralMath.mulDivF(FIXED_1, a, b)), FIXED_1);
}}
/**
* @dev Compute exp(a / b)
*/
function exp(uint256 a, uint256 b) internal pure returns (uint256, uint256) { unchecked {
return (fixedExp(IntegralMath.mulDivF(FIXED_1, a, b)), FIXED_1);
}}
/**
* @dev Compute log(x / FIXED_1) * FIXED_1
* Input range: FIXED_1 <= x <= 2 ^ 256 - 1
* Detailed description:
* - For x < LOG_MID, compute log(x)
* - For any other x, compute log(x / 2 ^ log2(x)) + log2(x) * log(2)
* - The value of log(2) is represented as floor(log(2) * FIXED_1)
* - With k = log2(x), this solution relies on the following identity:
* log(x) =
* log(x) + log(2 ^ k) - log(2 ^ k) =
* log(x) - log(2 ^ k) + log(2 ^ k) =
* log(x / log(2 ^ k)) + log(2 ^ k) =
* log(x / log(2 ^ k)) + k * log(2)
*/
function fixedLog(uint256 x) internal pure returns (uint256) { unchecked {
if (x < FIXED_1)
revert("fixedLog: x < min");
if (x < LOG_MID)
return optimalLog(x);
uint8 count = IntegralMath.floorLog2(x / FIXED_1);
return optimalLog(x >> count) + count * LN2_MIN;
}}
/**
* @dev Compute exp(x / FIXED_1) * FIXED_1
* Input range: 0 <= x <= EXP_MAX - 1
* Detailed description:
* - For x < EXP_MID, compute exp(x)
* - For any other x, compute exp(x % log(2)) * 2 ^ (x / log(2))
* - The value of log(2) is represented as ceil(log(2) * FIXED_1)
* - With k = x / log(2), this solution relies on the following identity:
* exp(x) =
* exp(x) * 2 ^ k / 2 ^ k =
* exp(x) * 2 ^ k / exp(k * log(2)) =
* exp(x) / exp(k * log(2)) * 2 ^ k =
* exp(x - k * log(2)) * 2 ^ k
*/
function fixedExp(uint256 x) internal pure returns (uint256) { unchecked {
if (x < EXP_MID)
return optimalExp(x);
if (x < EXP_MAX)
return optimalExp(x % LN2_MAX) << (x / LN2_MAX);
revert("fixedExp: x > max");
}}
/**
* @dev Compute log(x / FIXED_1) * FIXED_1
* Input range: FIXED_1 <= x <= FIXED_1 * 4 - 1
* Auto-generated via 'PrintAnalyticMathOptimalLog.py'
* Detailed description:
* - Rewrite the input as a product of natural exponents and a single residual r, such that 1 < r < 2
* - The natural logarithm of each (pre-calculated) exponent is the degree of the exponent
* - The natural logarithm of r is calculated via Taylor series for log(1 + x), where x = r - 1
* - The natural logarithm of the input is calculated by summing up the intermediate results above
* - For example: log(250) = log(e^4 * e^1 * e^0.5 * 1.021692859) = 4 + 1 + 0.5 + log(1 + 0.021692859)
*/
function optimalLog(uint256 x) internal pure returns (uint256) { unchecked {
uint256 res = 0;
uint256 y;
uint256 z;
uint256 w;
if (x > 0xd3094c70f034de4b96ff7d5b6f99fcd8) {res |= 0x40000000000000000000000000000000; x = x * 0x0bc5ab1b16779be3575bd8f0520a9d5b9 / 0x1368b2fc6f9609fe7aceb46aa619b8003;} // add 2^(-1)
if (x > 0xa45af1e1f40c333b3de1db4dd55f29a7) {res |= 0x20000000000000000000000000000000; x = x * 0x1368b2fc6f9609fe7aceb46aa619a2ef6 / 0x18ebef9eac820ae8682b9793ac6cffa93;} // add 2^(-2)
if (x > 0x910b022db7ae67ce76b441c27035c6a1) {res |= 0x10000000000000000000000000000000; x = x * 0x18ebef9eac820ae8682b9793ac6d11622 / 0x1c3d6a24ed82218787d624d3e5eb9a8c6;} // add 2^(-3)
if (x > 0x88415abbe9a76bead8d00cf112e4d4a8) {res |= 0x08000000000000000000000000000000; x = x * 0x1c3d6a24ed82218787d624d3e5eba8beb / 0x1e0fabfbc702a3ce5e31fe0609358b04d;} // add 2^(-4)
if (x > 0x84102b00893f64c705e841d5d4064bd3) {res |= 0x04000000000000000000000000000000; x = x * 0x1e0fabfbc702a3ce5e31fe06093589585 / 0x1f03f56a88b5d7914b00abf9776270f29;} // add 2^(-5)
if (x > 0x8204055aaef1c8bd5c3259f4822735a2) {res |= 0x02000000000000000000000000000000; x = x * 0x1f03f56a88b5d7914b00abf9776267973 / 0x1f80feabfeefa4927d10bdd54ead4eaf0;} // add 2^(-6)
if (x > 0x810100ab00222d861931c15e39b44e99) {res |= 0x01000000000000000000000000000000; x = x * 0x1f80feabfeefa4927d10bdd54ead599f9 / 0x1fc03fd56aa224f97fcbf133298883271;} // add 2^(-7)
if (x > 0x808040155aabbbe9451521693554f733) {res |= 0x00800000000000000000000000000000; x = x * 0x1fc03fd56aa224f97fcbf13329886bd10 / 0x1fe00ffaabffbbc71ad1e1184afc01529;} // add 2^(-8)
z = y = x - FIXED_1;
w = y * y / FIXED_1;
res += z * (0x100000000000000000000000000000000 - y) / 0x100000000000000000000000000000000; z = z * w / FIXED_1; // add y^01 / 01 - y^02 / 02
res += z * (0x0aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa - y) / 0x200000000000000000000000000000000; z = z * w / FIXED_1; // add y^03 / 03 - y^04 / 04
res += z * (0x099999999999999999999999999999999 - y) / 0x300000000000000000000000000000000; z = z * w / FIXED_1; // add y^05 / 05 - y^06 / 06
res += z * (0x092492492492492492492492492492492 - y) / 0x400000000000000000000000000000000; z = z * w / FIXED_1; // add y^07 / 07 - y^08 / 08
res += z * (0x08e38e38e38e38e38e38e38e38e38e38e - y) / 0x500000000000000000000000000000000; z = z * w / FIXED_1; // add y^09 / 09 - y^10 / 10
res += z * (0x08ba2e8ba2e8ba2e8ba2e8ba2e8ba2e8b - y) / 0x600000000000000000000000000000000; z = z * w / FIXED_1; // add y^11 / 11 - y^12 / 12
res += z * (0x089d89d89d89d89d89d89d89d89d89d89 - y) / 0x700000000000000000000000000000000; z = z * w / FIXED_1; // add y^13 / 13 - y^14 / 14
res += z * (0x088888888888888888888888888888888 - y) / 0x800000000000000000000000000000000; // add y^15 / 15 - y^16 / 16
return res;
}}
/**
* @dev Compute exp(x / FIXED_1) * FIXED_1
* Input range: 0 <= x <= EXP_MID - 1
* Auto-generated via 'PrintAnalyticMathOptimalExp.py'
* Detailed description:
* - Rewrite the input as a sum of binary exponents and a single residual r, as small as possible
* - The exponentiation of each binary exponent is given (pre-calculated)
* - The exponentiation of r is calculated via Taylor series for e^x, where x = r
* - The exponentiation of the input is calculated by multiplying the intermediate results above
* - For example: e^5.521692859 = e^(4 + 1 + 0.5 + 0.021692859) = e^4 * e^1 * e^0.5 * e^0.021692859
*/
function optimalExp(uint256 x) internal pure returns (uint256) { unchecked {
uint256 res = 0;
uint256 y;
uint256 z;
z = y = x % 0x10000000000000000000000000000000; // get the input modulo 2^(-3)
z = z * y / FIXED_1; res += z * 0x10e1b3be415a0000; // add y^02 * (20! / 02!)
z = z * y / FIXED_1; res += z * 0x05a0913f6b1e0000; // add y^03 * (20! / 03!)
z = z * y / FIXED_1; res += z * 0x0168244fdac78000; // add y^04 * (20! / 04!)
z = z * y / FIXED_1; res += z * 0x004807432bc18000; // add y^05 * (20! / 05!)
z = z * y / FIXED_1; res += z * 0x000c0135dca04000; // add y^06 * (20! / 06!)
z = z * y / FIXED_1; res += z * 0x0001b707b1cdc000; // add y^07 * (20! / 07!)
z = z * y / FIXED_1; res += z * 0x000036e0f639b800; // add y^08 * (20! / 08!)
z = z * y / FIXED_1; res += z * 0x00000618fee9f800; // add y^09 * (20! / 09!)
z = z * y / FIXED_1; res += z * 0x0000009c197dcc00; // add y^10 * (20! / 10!)
z = z * y / FIXED_1; res += z * 0x0000000e30dce400; // add y^11 * (20! / 11!)
z = z * y / FIXED_1; res += z * 0x000000012ebd1300; // add y^12 * (20! / 12!)
z = z * y / FIXED_1; res += z * 0x0000000017499f00; // add y^13 * (20! / 13!)
z = z * y / FIXED_1; res += z * 0x0000000001a9d480; // add y^14 * (20! / 14!)
z = z * y / FIXED_1; res += z * 0x00000000001c6380; // add y^15 * (20! / 15!)
z = z * y / FIXED_1; res += z * 0x000000000001c638; // add y^16 * (20! / 16!)
z = z * y / FIXED_1; res += z * 0x0000000000001ab8; // add y^17 * (20! / 17!)
z = z * y / FIXED_1; res += z * 0x000000000000017c; // add y^18 * (20! / 18!)
z = z * y / FIXED_1; res += z * 0x0000000000000014; // add y^19 * (20! / 19!)
z = z * y / FIXED_1; res += z * 0x0000000000000001; // add y^20 * (20! / 20!)
res = res / 0x21c3677c82b40000 + y + FIXED_1; // divide by 20! and then add y^1 / 1! + y^0 / 0!
if ((x & 0x010000000000000000000000000000000) != 0) res = res * 0x1c3d6a24ed82218787d624d3e5eba2a0f / 0x18ebef9eac820ae8682b9793ac6d1883a; // multiply by e^2^(-3)
if ((x & 0x020000000000000000000000000000000) != 0) res = res * 0x18ebef9eac820ae8682b9793ac6d1e726 / 0x1368b2fc6f9609fe7aceb46aa619bae94; // multiply by e^2^(-2)
if ((x & 0x040000000000000000000000000000000) != 0) res = res * 0x1368b2fc6f9609fe7aceb46aa6199a7d7 / 0x0bc5ab1b16779be3575bd8f0520a8b756; // multiply by e^2^(-1)
if ((x & 0x080000000000000000000000000000000) != 0) res = res * 0x0bc5ab1b16779be3575bd8f0520a67cd2 / 0x0454aaa8efe072e7f6ddbab84b409101a; // multiply by e^2^(+0)
if ((x & 0x100000000000000000000000000000000) != 0) res = res * 0x0454aaa8efe072e7f6ddbab84b408736f / 0x00960aadc109e7a3bf45780996156d0a3; // multiply by e^2^(+1)
if ((x & 0x200000000000000000000000000000000) != 0) res = res * 0x00960aadc109e7a3bf45780996149ade3 / 0x0002bf84208204f5977f9a8cf01fd8f74; // multiply by e^2^(+2)
return res;
}}
}