|
| 1 | +// SPDX-License-Identifier: MIT |
| 2 | +pragma solidity ^0.8.24; |
| 3 | + |
| 4 | +import {SD59x18, convert, ZERO} from "@prb/math/SD59x18.sol"; |
| 5 | +import {SQRT_PI, SQRT_TWO_PI, EPSILON, DUST_THRESHOLD, NORM_TOLERANCE} from "../types/Types.sol"; |
| 6 | +import {InvalidNorm, NotCriticalPoint, NotMinimum, WrongSide, DustPosition} from "../types/Types.sol"; |
| 7 | + |
| 8 | +/// @title GaussianMath |
| 9 | +/// @notice Pure math library for Gaussian PDF evaluation and verification. |
| 10 | +/// @dev All functions use PRBMath SD59x18 fixed-point arithmetic. |
| 11 | +/// The Gaussian is f(x) = λ/(σ·√(2π)) · exp(−(x−μ)²/(2σ²)) |
| 12 | +/// i.e. f = λ · φ(x; μ, σ) where φ is the standard Normal PDF. |
| 13 | +library GaussianMath { |
| 14 | + /// @notice Evaluate f(x) = λ/(σ·√(2π)) · exp(−(x−μ)²/(2σ²)) |
| 15 | + function eval(SD59x18 mu, SD59x18 lambda, SD59x18 sigma, SD59x18 x) internal pure returns (SD59x18) { |
| 16 | + SD59x18 amplitude = lambda / (sigma * SQRT_TWO_PI); |
| 17 | + return amplitude * _expTerm(mu, sigma, x); |
| 18 | + } |
| 19 | + |
| 20 | + /// @notice First derivative: f'(x) = −(x−μ)/σ² · f(x) |
| 21 | + function derivative(SD59x18 mu, SD59x18 lambda, SD59x18 sigma, SD59x18 x) internal pure returns (SD59x18) { |
| 22 | + SD59x18 diff = mu - x; |
| 23 | + SD59x18 sigmaSq = sigma * sigma; |
| 24 | + return (diff / sigmaSq) * eval(mu, lambda, sigma, x); |
| 25 | + } |
| 26 | + |
| 27 | + /// @notice Second derivative: f''(x) = ((x−μ)²/σ⁴ − 1/σ²) · f(x) |
| 28 | + function secondDerivative(SD59x18 mu, SD59x18 lambda, SD59x18 sigma, SD59x18 x) |
| 29 | + internal |
| 30 | + pure |
| 31 | + returns (SD59x18) |
| 32 | + { |
| 33 | + SD59x18 diff = x - mu; |
| 34 | + SD59x18 sigmaSq = sigma * sigma; |
| 35 | + SD59x18 factor = (diff * diff) / (sigmaSq * sigmaSq) - convert(1) / sigmaSq; |
| 36 | + return factor * eval(mu, lambda, sigma, x); |
| 37 | + } |
| 38 | + |
| 39 | + /// @notice Verify the L₂ norm invariant: λ² ≈ 2·k²·σ·√π |
| 40 | + /// @dev Uses relative tolerance NORM_TOLERANCE (0.1%). |
| 41 | + function verifyNorm(SD59x18 lambda, SD59x18 sigma, SD59x18 kSquared) internal pure { |
| 42 | + SD59x18 lhs = lambda * lambda; |
| 43 | + SD59x18 rhs = convert(2) * kSquared * sigma * SQRT_PI; |
| 44 | + |
| 45 | + SD59x18 diff = lhs > rhs ? lhs - rhs : rhs - lhs; |
| 46 | + SD59x18 maxVal = lhs > rhs ? lhs : rhs; |
| 47 | + |
| 48 | + if (diff / maxVal > NORM_TOLERANCE) revert InvalidNorm(); |
| 49 | + } |
| 50 | + |
| 51 | + /// @notice 4-check min-point verification for collateral calculation. |
| 52 | + /// @param fMu Current aggregate mean |
| 53 | + /// @param fLambda Current aggregate lambda |
| 54 | + /// @param fSigma Current aggregate sigma |
| 55 | + /// @param gMu Proposed new aggregate mean |
| 56 | + /// @param gLambda Proposed new aggregate lambda |
| 57 | + /// @param gSigma Proposed new aggregate sigma |
| 58 | + /// @param candidate The x-coordinate of the claimed minimum of g(x)−f(x) |
| 59 | + /// @return minValue The (negative) value g(candidate)−f(candidate); collateral = −minValue |
| 60 | + function verifyMinPoint( |
| 61 | + SD59x18 fMu, |
| 62 | + SD59x18 fLambda, |
| 63 | + SD59x18 fSigma, |
| 64 | + SD59x18 gMu, |
| 65 | + SD59x18 gLambda, |
| 66 | + SD59x18 gSigma, |
| 67 | + SD59x18 candidate |
| 68 | + ) internal pure returns (SD59x18 minValue) { |
| 69 | + // Check 1: Critical point — |g'(x) − f'(x)| < EPSILON |
| 70 | + SD59x18 gPrime = derivative(gMu, gLambda, gSigma, candidate); |
| 71 | + SD59x18 fPrime = derivative(fMu, fLambda, fSigma, candidate); |
| 72 | + SD59x18 firstDerivDiff = gPrime - fPrime; |
| 73 | + if (firstDerivDiff < ZERO) firstDerivDiff = -firstDerivDiff; |
| 74 | + if (firstDerivDiff > EPSILON) revert NotCriticalPoint(); |
| 75 | + |
| 76 | + // Check 2: Minimum — g''(x) − f''(x) > 0 |
| 77 | + SD59x18 gDoublePrime = secondDerivative(gMu, gLambda, gSigma, candidate); |
| 78 | + SD59x18 fDoublePrime = secondDerivative(fMu, fLambda, fSigma, candidate); |
| 79 | + if (gDoublePrime - fDoublePrime < ZERO) revert NotMinimum(); |
| 80 | + |
| 81 | + // Check 3: Geometric side check — candidate must be on opposite side of g's mean from f's mean |
| 82 | + if (gMu > fMu) { |
| 83 | + if (candidate >= fMu) revert WrongSide(); |
| 84 | + } else if (gMu < fMu) { |
| 85 | + if (candidate <= fMu) revert WrongSide(); |
| 86 | + } else { |
| 87 | + // Equal means: minimum is in the tails, candidate must be at least 1σ from shared mean |
| 88 | + SD59x18 dist = candidate - fMu; |
| 89 | + if (dist < ZERO) dist = -dist; |
| 90 | + if (dist < fSigma) revert WrongSide(); |
| 91 | + } |
| 92 | + |
| 93 | + // Value computation (for Check 4 and return value) |
| 94 | + SD59x18 gVal = eval(gMu, gLambda, gSigma, candidate); |
| 95 | + SD59x18 fVal = eval(fMu, fLambda, fSigma, candidate); |
| 96 | + minValue = gVal - fVal; |
| 97 | + |
| 98 | + // Check 4: Dust threshold — min value is meaningfully negative |
| 99 | + if (minValue > -DUST_THRESHOLD) revert DustPosition(); |
| 100 | + } |
| 101 | + |
| 102 | + /// @notice Derive lambda from the invariant: λ = √(2·k²·σ·√π) |
| 103 | + function computeLambda(SD59x18 kSquared, SD59x18 sigma) internal pure returns (SD59x18) { |
| 104 | + SD59x18 inner = convert(2) * kSquared * sigma * SQRT_PI; |
| 105 | + return inner.sqrt(); |
| 106 | + } |
| 107 | + |
| 108 | + /// @notice Closed-form min of g(x)−f(x) for cross-checking in tests. |
| 109 | + /// @dev m = λ_g/(2·σ_g·√π) − λ_f/√(2π·(σ_f²+σ_g²)) · exp(−(μ_f−μ_g)²/(2·(σ_f²+σ_g²))) |
| 110 | + function computeFormulaMin( |
| 111 | + SD59x18 fMu, |
| 112 | + SD59x18 fLambda, |
| 113 | + SD59x18 fSigma, |
| 114 | + SD59x18 gMu, |
| 115 | + SD59x18 gLambda, |
| 116 | + SD59x18 gSigma |
| 117 | + ) internal pure returns (SD59x18) { |
| 118 | + // Term 1: Self-norm — λ_g / (2·σ_g·√π) |
| 119 | + SD59x18 selfNorm = gLambda / (convert(2) * gSigma * SQRT_PI); |
| 120 | + |
| 121 | + // Combined sigma for the cross-term |
| 122 | + SD59x18 fSigmaSq = fSigma * fSigma; |
| 123 | + SD59x18 gSigmaSq = gSigma * gSigma; |
| 124 | + SD59x18 combinedSigmaSq = fSigmaSq + gSigmaSq; |
| 125 | + SD59x18 combinedSigma = combinedSigmaSq.sqrt(); |
| 126 | + |
| 127 | + // Cross-term amplitude: λ_f / (√(2π) · combinedSigma) |
| 128 | + SD59x18 crossAmplitude = fLambda / (SQRT_TWO_PI * combinedSigma); |
| 129 | + |
| 130 | + // Cross-term exponent: −(μ_f − μ_g)² / (2·(σ_f² + σ_g²)) |
| 131 | + SD59x18 muDiff = fMu - gMu; |
| 132 | + SD59x18 exponent = -(muDiff * muDiff) / (convert(2) * combinedSigmaSq); |
| 133 | + SD59x18 crossTerm = crossAmplitude * exponent.exp(); |
| 134 | + |
| 135 | + return selfNorm - crossTerm; |
| 136 | + } |
| 137 | + |
| 138 | + /// @dev Shared exponential term: exp(−(x−μ)²/(2σ²)) |
| 139 | + function _expTerm(SD59x18 mu, SD59x18 sigma, SD59x18 x) private pure returns (SD59x18) { |
| 140 | + SD59x18 diff = x - mu; |
| 141 | + SD59x18 exponent = -(diff * diff) / (convert(2) * sigma * sigma); |
| 142 | + return exponent.exp(); |
| 143 | + } |
| 144 | +} |
0 commit comments