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