72 lines
2.5 KiB
Solidity
72 lines
2.5 KiB
Solidity
import './ABDKMathQuad.sol';
|
|
|
|
|
|
contract pi {
|
|
using ABDKMathQuad for bytes16;
|
|
|
|
// Immutable variables (set once in constructor)
|
|
bytes16 public immutable C_426880;
|
|
bytes16 public immutable C_10005;
|
|
bytes16 public immutable C_13591409;
|
|
bytes16 public immutable C_545140134;
|
|
bytes16 public immutable C_640320;
|
|
bytes16 public immutable C_12;
|
|
|
|
constructor() {
|
|
C_426880 = ABDKMathQuad.fromInt(426880);
|
|
C_10005 = ABDKMathQuad.fromInt(10005);
|
|
C_13591409 = ABDKMathQuad.fromUInt(13591409);
|
|
C_545140134 = ABDKMathQuad.fromUInt(545140134);
|
|
C_640320 = ABDKMathQuad.fromUInt(640320);
|
|
C_12 = ABDKMathQuad.fromUInt(12);
|
|
}
|
|
|
|
// Compute factorial of n (as uint256), note: limited by gas
|
|
function factorial(uint256 n) internal pure returns (uint256) {
|
|
if (n == 0 || n == 1) return 1;
|
|
uint256 result = 1;
|
|
for (uint256 i = 2; i <= n; i++) {
|
|
result *= i;
|
|
}
|
|
return result;
|
|
}
|
|
|
|
// Compute power (uint256 base ^ uint256 exp)
|
|
function pow(uint256 base, uint256 exp) internal pure returns (uint256) {
|
|
uint256 result = 1;
|
|
for (uint256 i = 0; i < exp; i++) {
|
|
result *= base;
|
|
}
|
|
return result;
|
|
}
|
|
|
|
// Compute one term of the Chudnovsky series for k
|
|
function chudnovskyTerm(uint256 k) internal pure returns (bytes16 numerator, bytes16 denominator) {
|
|
uint256 sixKFact = factorial(6 * k);
|
|
uint256 kFact = factorial(k);
|
|
uint256 threeKFact = factorial(3 * k);
|
|
|
|
// Use int256 to allow negative multiplication
|
|
int256 numeratorInt = int256(sixKFact) * int256(13591409 + 545140134 * k);
|
|
if (k % 2 == 1) numeratorInt *= -1; // Correctly applies sign
|
|
|
|
uint256 denominatorInt = threeKFact * (kFact ** 3) * pow(640320, 3 * k);
|
|
|
|
numerator = ABDKMathQuad.fromInt(numeratorInt); // Ensure ABDKMathQuad supports int
|
|
denominator = ABDKMathQuad.fromUInt(denominatorInt);
|
|
}
|
|
|
|
// Approximate pi using n terms (WARNING: only small n due to gas and uint256 limits)
|
|
function computePi(uint256 n) public view returns (bytes16) {
|
|
bytes16 sum = ABDKMathQuad.fromUInt(0);
|
|
|
|
for (uint256 k = 0; k < n; k++) {
|
|
(bytes16 num, bytes16 den) = chudnovskyTerm(k);
|
|
sum = sum.add(num.div(den));
|
|
}
|
|
|
|
bytes16 sqrt10005 = ABDKMathQuad.sqrt(C_10005);
|
|
bytes16 factor = C_426880.mul(sqrt10005);
|
|
return factor.div(sum);
|
|
}
|
|
} |