fix: input range of current Math.cdf implementation

[dlspn]

Summary

With the current implementation of cdf and ncdf, the input range with correct output is a bit small due to overflows and the algorithm to calculate exponentials in optimalExp. This pull request fixes the problems and the modified cdf will return good results for all possible inputs with reasonable precision.

Problem

rvol/contracts/libraries/Math.sol

Line 247 in ee066ee

uint256 exp = ((x / 2) * x) / FIXED_1;

In Math.sol L247 and L267, the calculation of variable exp overflows when $$\log_2{\frac{x^2}{2}}\geq256$$ i.e. $$x\geq2^{128.5}$$ considering the variable x is 127-bit binary fixed-point, the actual range is $$x\geq2^{128.5}\div2^{127}=2^{1.5}\approx2.8284$$
which means the function can only return correct results for x < 2.82, approximately:

In rvol we normally pass the Black-Scholes parameters d1 and d2 to the CDF functions. They are calculated from spot price, strike, volatility and time to expire of the option or underlying asset. With a reasonable set of these parameters

• spot price is 1.65x of strike
• annualized volatility is 108%
• 7 days left to expire

we get d1 = 3.408 approximately, which is far larger than the 2.82 limit. So cdf can give wrong answers sometimes in options-related contexts.

Solution

First we replace the uint256 exp = ((x / 2) * x) / FIXED_1; line
with uint256 exp = (x * (x >> 3)) / (FIXED_1 >> 2);
The 3-bit shift is enough to prevent overflow for $$x\lt\sqrt{32}\approx5.6568542$$ (will explain it later)

After this modification the output of cdf is like

Now the input range extends to ~5 and it's enough for real options calculations. But we go a step further to make the CDF functions better.

rvol/contracts/libraries/Math.sol

Line 181 in ee066ee

if ((x & 0x400000000000000000000000000000000) != 0)

This time the periodic falls are due to the algorithm to implement optimalExp. After get Taylor series of small residuals it elaborates the exponential for x bit by bit but only to the 131st bit, which corresbonds to the 2^3 bit - so optimalExp will return periodic results with 16 as the period (in fixed-point).

Since the argument we pass to optimalExp is $\frac{x^2}{2}$, the output falls at $$\frac{x^2}{2}=16\cdot n\quad(n=1,2,3,\ldots)$$
thus the first falling point $$x_1=\sqrt{32}$$
the actual CDF value (0.9999999...) has been very close to 1 since long before this point so we let the functions return 1 if input is greater than it (5.656854)

if (x > (5656854 * FIXED_1) / 1e6) {
    return 1e14;
}

from previous discussions we know the input exp will not overflow either within this range so cdf can return good results now:

One more thing

rvol/contracts/libraries/Math.sol

Line 83 in ee066ee

if (xAux >= 0x8) {

the constant in L83 in sqrt should be 0x4 for the last bit in the guess of initial value.