From 54b6292077af418c903dce6efa7cb2d0a20d8ea2 Mon Sep 17 00:00:00 2001
From: alex <alexangel.j@gmail.com>
Date: Sat, 20 May 2023 22:33:42 -0700
Subject: [PATCH 1/9] feat(c-port): normal cdf with 1e18 precision

---
 src/Ndtr.sol                    | 219 +++++++++++++++++++++++++++++++
 src/Ndtri.sol                   | 187 +++++++++++++++++++++++++++
 src/Pdtr.sol                    |  11 ++
 src/libraries/ArrayLib.sol      |  56 ++++++++
 src/libraries/PolynomialLib.sol |  49 +++++++
 src/libraries/RayMathLib.sol    |  64 ++++++++++
 src/test/Ndtr.t.sol             | 220 ++++++++++++++++++++++++++++++++
 src/test/Ndtri.t.sol            | 145 +++++++++++++++++++++
 src/test/RayMathLib.t.sol       |  15 +++
 9 files changed, 966 insertions(+)
 create mode 100644 src/Ndtr.sol
 create mode 100644 src/Ndtri.sol
 create mode 100644 src/Pdtr.sol
 create mode 100644 src/libraries/ArrayLib.sol
 create mode 100644 src/libraries/PolynomialLib.sol
 create mode 100644 src/libraries/RayMathLib.sol
 create mode 100644 src/test/Ndtr.t.sol
 create mode 100644 src/test/Ndtri.t.sol
 create mode 100644 src/test/RayMathLib.t.sol

diff --git a/src/Ndtr.sol b/src/Ndtr.sol
new file mode 100644
index 0000000..67c076f
--- /dev/null
+++ b/src/Ndtr.sol
@@ -0,0 +1,219 @@
+// SPDX-License-Identifier: AGPL-3.0-only
+pragma solidity ^0.8.4;
+
+import "solmate/utils/FixedPointMathLib.sol";
+
+import "./libraries/RayMathLib.sol";
+import "./libraries/PolynomialLib.sol";
+import "./libraries/ArrayLib.sol";
+
+import {console2 as logger} from "forge-std/console2.sol"; // todo: remove
+
+/// @dev This library implements the Normal cumulative distribution function (CDF).
+///      The implementation is inspired by the algorithm implemented in the Cephes library.
+/// @custom:reference https://github.com/jeremybarnes/cephes/blob/master/cprob/ndtr.c
+library Ndtr {
+    ////////////////
+    // Constants  //
+    ////////////////
+
+    /// @dev
+    /// In the context of the erf and erfc functions,
+    /// the square root of 0.5 is used as a boundary value for a range check.
+    /// When the absolute value of the input x is smaller than the square root of 0.5,
+    /// a series expansion is used to compute the function's value.
+    /// Otherwise, a continued fraction expansion is used.
+    ///
+    /// These two different methods for computing the error function and
+    /// the complementary error function provide a balance between
+    /// accuracy and computational efficiency,
+    /// depending on the size of the input.
+    int256 constant RAY_SQRTH = 0.7071067811865475244e27; // sqrt(5) = 0.7071067811865475244
+
+    uint256 constant P_LENGTH = 9;
+    bytes constant P = abi.encodePacked(
+        [
+            int256(2.46196981473530512e17), // P[0] is missing 524, decimal places 28-30. todo: check out?
+            int256(5.64189564831068821977e26),
+            int256(7.46321056442269912687e27),
+            int256(4.86371970985681366614e28),
+            int256(1.96520832956077098242e29),
+            int256(5.26445194995477358631e29),
+            int256(9.3452852717195760754e29),
+            int256(1.02755188689515710272e30),
+            int256(5.57535335369399327526e29)
+        ]
+    );
+
+    uint256 constant Q_LENGTH = 8;
+    bytes constant Q = abi.encodePacked(
+        [
+            // 1.0e27
+            1.32281951154744992508e28,
+            8.67072140885989742329e28,
+            3.54937778887819891062e29,
+            9.75708501743205489753e29,
+            1.82390916687909736289e30,
+            2.24633760818710981792e30,
+            1.65666309194161350182e30,
+            5.57535340817727675546e29
+        ]
+    );
+
+    uint256 constant R_LENGTH = 6;
+    bytes constant R = abi.encodePacked(
+        [
+            5.64189583547755073984e26,
+            1.27536670759978104416e27,
+            5.01905042251180477414e27,
+            6.16021097993053585195e27,
+            7.4097426995044893916e27,
+            2.9788666537210024067e27
+        ]
+    );
+
+    uint256 constant S_LENGTH = 6;
+    bytes constant S = abi.encodePacked(
+        [
+            // 1.0e27
+            2.2605286322011727659e27,
+            9.39603524938001434673e27,
+            1.20489539808096656605e28,
+            1.70814450747565897222e28,
+            9.60896809063285878198e27,
+            3.3690764510008151605e27
+        ]
+    );
+
+    uint256 constant T_LENGTH = 5;
+    bytes constant T = abi.encodePacked(
+        [
+            9.60497373987051638749e27,
+            9.00260197203842689217e28,
+            2.23200534594684319226e30,
+            7.00332514112805075473e30,
+            5.55923013010394962768e31
+        ]
+    );
+
+    uint256 constant U_LENGTH = 5;
+    bytes constant U = abi.encodePacked(
+        [
+            // 1.0e27
+            3.35617141647503099647e28,
+            5.21357949780152679795e29,
+            4.59432382970980127987e30,
+            2.26290000613890934246e31,
+            4.92673942608635921086e31
+        ]
+    );
+
+    ////////////////
+    // Functions  //
+    ////////////////
+
+    /// todo: need to handling rounding, unless we want to default to truncation.
+    function ndtr(int256 a) internal view returns (int256) {
+        int256 x;
+        int256 y;
+        int256 z;
+
+        x = mulfp(a, RAY_SQRTH); // x = a * RAY_SQRTH;
+        z = absolute(x);
+
+        if (z < RAY_ONE) {
+            y = RAY_HALF + mulfp(RAY_HALF, erf(x)); // y = RAY_HALF + RAY_HALF * erf(x);
+        } else {
+            y = mulfp(RAY_HALF, erfc(z)); // y = RAY_HALF * erfc(z);
+
+            if (x > 0) y = RAY_ONE - y;
+        }
+
+        return y;
+    }
+
+    function erfc(int256 a) internal view returns (int256) {
+        int256 p;
+        int256 q;
+        int256 x;
+        int256 y;
+        int256 z;
+
+        if (a < 0) x = -a;
+        else x = a;
+
+        if (x < RAY_ONE) return RAY_ONE - erf(a);
+
+        z = -mulfp(a, a); // z = -a * a;
+
+        if (z < -RAY_MAXLOG) {
+            if (a < 0) return RAY_TWO;
+            else return 0;
+        }
+
+        z = FixedPointMathLib.expWad(z / 1e9) * 1e9;
+
+        if (x < RAY_EIGHT) {
+            int256[] memory input0 = copy9(abi.decode(P, (int256[9])));
+            p = polevl(x, input0, 8); // p = P[0] + P[1]*x + ... + P[8]*x^8
+            int256[] memory input1 = copy8(abi.decode(Q, (int256[8])));
+            q = p1evl(x, input1, 8); // q = 1 + Q[0]*x + ... + Q[7]*x^7
+        } else {
+            int256[] memory input0 = copy6(abi.decode(R, (int256[6])));
+            p = polevl(x, input0, 5); // p = R[0] + R[1]*x + ... + R[5]*x^5
+            int256[] memory input1 = copy6(abi.decode(S, (int256[6])));
+            q = p1evl(x, input1, 6); // q = 1 + S[0]*x + ... + S[6]*x^6
+        }
+
+        y = (z * p) / q;
+
+        if (a < 0) y = RAY_TWO - y;
+
+        if (y == 0) {
+            if (a < 0) y = RAY_TWO;
+            else y = 0;
+        }
+
+        return y;
+    }
+
+    /// @dev Exponentially scaled erfc function.
+    ///      exp(x^2) erfc(x)
+    ///      valid for x > 1
+    function erfce(int256 x) internal view returns (int256) {
+        int256 p;
+        int256 q;
+
+        if (x < RAY_EIGHT) {
+            int256[] memory input0 = copy9(abi.decode(P, (int256[9])));
+            p = polevl(x, input0, 8); // p = P[0] + P[1]*x + ... + P[8]*x^8
+            int256[] memory input1 = copy8(abi.decode(Q, (int256[8])));
+            q = p1evl(x, input1, 8); // q = 1 + Q[0]*x + ... + Q[7]*x^7
+        } else {
+            int256[] memory input0 = copy6(abi.decode(R, (int256[6])));
+            p = polevl(x, input0, 5); // p = R[0] + R[1]*x + ... + R[5]*x^5
+            int256[] memory input1 = copy6(abi.decode(S, (int256[6])));
+            q = p1evl(x, input1, 6); // q = 1 + S[0]*x + ... + S[6]*x^6
+        }
+
+        return divfp(p, q); // p / q;
+    }
+
+    function erf(int256 x) internal view returns (int256) {
+        int256 y;
+        int256 z;
+
+        if (absolute(x) > RAY_ONE) return RAY_ONE - erfc(x);
+
+        z = mulfp(x, x); // z = x * x;
+
+        int256[] memory input0 = copy5(abi.decode(T, (int256[5])));
+        int256[] memory input1 = copy5(abi.decode(U, (int256[5])));
+        // p = T[0] + T[1]*z + ... + T[4]*z^4
+        // q = 1 + U[0]*z + ... + U[4]*z^5
+        // y = z * p / q
+        y = x * polevl(z, input0, 4) / p1evl(z, input1, 5);
+
+        return y;
+    }
+}
diff --git a/src/Ndtri.sol b/src/Ndtri.sol
new file mode 100644
index 0000000..ceb73e3
--- /dev/null
+++ b/src/Ndtri.sol
@@ -0,0 +1,187 @@
+// SPDX-License-Identifier: AGPL-3.0-only
+pragma solidity ^0.8.4;
+
+import "solmate/utils/FixedPointMathLib.sol";
+
+import "./libraries/RayMathLib.sol";
+import "./libraries/PolynomialLib.sol";
+import "./libraries/ArrayLib.sol";
+
+import {console2 as logger} from "forge-std/console2.sol"; // todo: remove
+
+/// @dev This library implements the Inverse normal cumulative distribution function (PPF).
+///      Also known as the percent point function or quantile function.
+///      The implementation is inspired by the algorithm implemented in the Cephes library.
+/// @custom:reference https://github.com/jeremybarnes/cephes/blob/master/cprob/ndtri.c
+library Ndtri {
+    ////////////////
+    // Errors     //
+    ////////////////
+
+    error MaxNumError();
+
+    //////////////////////////////////////////////
+    // Approximation for 0 <= |y - 0.5| <= 3/8  //
+    //////////////////////////////////////////////
+
+    int256 constant P0_LENGTH = 5;
+    bytes constant P0 = abi.encodePacked(
+        [
+            -int256(5.99633501014107895267e28),
+            int256(9.80010754185999661536e28),
+            -int256(5.66762857469070293439e28),
+            int256(1.39312609387279679503e28),
+            -int256(1.23916583867381258016e27)
+        ]
+    );
+
+    int256 constant Q0_LENGTH = 8;
+    bytes constant Q0 = abi.encodePacked(
+        [
+            int256(1.95448858338141759834e27),
+            int256(4.67627912898881538453e27),
+            int256(8.63602421390890590575e28),
+            -int256(2.25462687854119370527e29),
+            int256(2.00260212380060660359e29),
+            -int256(8.20372256168333339912e28),
+            int256(1.59056225126211695515e28),
+            -int256(1.18331621121330003142e27)
+        ]
+    );
+
+    /////////////////////////////////////////////////////////////////////
+    // Approximation for interval z = sqrt(-2 log y ) between 2 and 8  //
+    // i.e., y between exp(-2) = .135 and exp(-32) = 1.27e-14.         //
+    /////////////////////////////////////////////////////////////////////
+
+    uint256 constant P1_LENGTH = 9;
+    bytes constant P1 = abi.encodePacked(
+        [
+            int256(4.05544892305962419923e27),
+            int256(3.15251094599893866154e28),
+            int256(5.71628192246421288162e28),
+            int256(4.408050738932008347e28),
+            int256(1.46849561928858024014e28),
+            int256(2.18663306850790267539e27),
+            -int256(1.40256079171354495875e26),
+            -int256(3.50424626827848203418e25),
+            -int256(8.57456785154685413611e23)
+        ]
+    );
+
+    uint256 constant Q1_LENGTH = 8;
+    bytes constant Q1 = abi.encodePacked(
+        [
+            int256(1.57799883256466749731e28),
+            int256(4.53907635128879210584e28),
+            int256(4.1317203825467203044e28),
+            int256(1.50425385692907503408e28),
+            int256(2.50464946208309415979e27),
+            -int256(1.42182922854787788574e26),
+            -int256(3.80806407691578277194e25),
+            -int256(9.33259480895457427372e23)
+        ]
+    );
+
+    //////////////////////////////////////////////////////////////////////
+    // Approximation for interval z = sqrt(-2 log y ) between 8 and 64  //
+    // i.e., y between exp(-32) = 1.27e-14 and exp(-2048) = 3.67e-890.  //
+    //////////////////////////////////////////////////////////////////////
+
+    uint256 constant P2_LENGTH = 9;
+    bytes constant P2 = abi.encodePacked(
+        [
+            int256(3.2377489177694603597e27),
+            int256(6.91522889068984211695e27),
+            int256(3.93881025292474443415e27),
+            int256(1.33303460815807542389e27),
+            int256(2.01485389549179081538e26),
+            int256(1.23716634817820021358e25),
+            int256(3.01581553508235416007e23),
+            int256(2.65806974686737550832e21),
+            int256(6.239745391849832937e18) // todo: missing a 3 at the end! bad?
+        ]
+    );
+
+    uint256 constant Q2_LENGTH = 8;
+    bytes constant Q2 = abi.encodePacked(
+        [
+            int256(6.02427039364742014255e27),
+            int256(3.67983563856160859403e27),
+            int256(1.37702099489081330271e27),
+            int256(2.1623699359449663589e26),
+            int256(1.34204006088543189037e25),
+            int256(3.28014464682127739104e23),
+            int256(2.89247864745380683936e21),
+            int256(6.790194080099812744e18) // missing a 25 at the end! todo: fix?
+        ]
+    );
+
+    ////////////////
+    // Functions  //
+    ////////////////
+
+    /// todo: need sqrt function with enough precision
+    /// along with log
+    /// todo: precision. Looks like it loses precision towards the bounds. Maybe better sqrt resolves.
+    function ndtri(int256 y0) internal view returns (int256) {
+        int256 x;
+        int256 y;
+        int256 z;
+        int256 y2;
+        int256 x0;
+        int256 x1;
+
+        if (y0 <= 0) revert MaxNumError();
+        if (y0 >= RAY_ONE) revert MaxNumError();
+
+        int256 code = 1;
+
+        y = y0;
+        if (y > RAY_ONE - 0.13533528323661269189e27) {
+            // 0.135... = exp(-2)
+            y = RAY_ONE - y;
+            code = 0;
+        }
+
+        if (y > 0.13533528323661269189e27) {
+            int256[] memory P0_ARRAY = copy5(abi.decode(P0, (int256[5])));
+            int256[] memory Q0_ARRAY = copy8(abi.decode(Q0, (int256[8])));
+            y = y - RAY_HALF;
+            y2 = mulfp(y, y);
+            x = y + mulfp(y, (y2 * polevl(y2, P0_ARRAY, 4) / p1evl(y2, Q0_ARRAY, 8)));
+            x = mulfp(x, RAY_SQRT2PI);
+            return x;
+        }
+
+        // `y` is assumed to be a probability between 0 and 1.
+        // This transformation is part of the formula for the quantile function of the normal distribution
+        // and serves to adjust the scale and distribution of the data.
+        // It also helps ensure that x will be a positive number because we're taking the square root.
+        // The -2.0 * log(y) essentially gives the quantile of the exponential distribution,
+        // which is then square rooted to approximate the quantile for the normal distribution.
+        x = sqrtfp(mulfp(-RAY_TWO, logfp(y))); // √(-2 * ln(y))
+
+        // Part of a numerical approximation method (like a variant of Newton's method)
+        // to refine the value of x to make it a better approximation.
+        x0 = x - divfp(logfp(x), x); // x - ln(x) / x
+
+        z = divfp(RAY_ONE, x);
+        if (x < RAY_EIGHT) {
+            int256[] memory P1_ARRAY = copy9(abi.decode(P1, (int256[9])));
+            int256[] memory Q1_ARRAY = copy8(abi.decode(Q1, (int256[8])));
+            x1 = z * polevl(z, P1_ARRAY, uint256(8)) / p1evl(z, Q1_ARRAY, uint256(8));
+        } else {
+            int256[] memory P2_ARRAY = copy9(abi.decode(P2, (int256[9])));
+            int256[] memory Q2_ARRAY = copy8(abi.decode(Q2, (int256[8])));
+            x1 = z * polevl(z, P2_ARRAY, uint256(8)) / p1evl(z, Q2_ARRAY, uint256(8));
+        }
+
+        x = x0 - y;
+        if (code != 0) {
+            x = -x;
+        }
+
+        return x;
+    }
+}
diff --git a/src/Pdtr.sol b/src/Pdtr.sol
new file mode 100644
index 0000000..7f67ebe
--- /dev/null
+++ b/src/Pdtr.sol
@@ -0,0 +1,11 @@
+// SPDX-License-Identifier: AGPL-3.0-only
+pragma solidity ^0.8.4;
+
+import "solmate/utils/FixedPointMathLib.sol";
+
+import "./libraries/RayMathLib.sol";
+import "./libraries/PolynomialLib.sol";
+import "./libraries/ArrayLib.sol";
+
+/// todo
+library Pdtr {}
diff --git a/src/libraries/ArrayLib.sol b/src/libraries/ArrayLib.sol
new file mode 100644
index 0000000..eed01a9
--- /dev/null
+++ b/src/libraries/ArrayLib.sol
@@ -0,0 +1,56 @@
+// SPDX-License-Identifier: GPL-3.0-only
+pragma solidity ^0.8.4;
+
+////////////////
+// Functions  //
+////////////////
+
+/// @dev Copies a fixed-size array of 5 elements into a new dynamically allocated array.
+function copy5(int256[5] memory input) pure returns (int256[] memory) {
+    uint256 len = input.length;
+    int256[] memory output = new int256[](len);
+    for (uint256 i; i < len; i++) {
+        output[i] = input[i];
+    }
+    return output;
+}
+
+/// @dev Copies a fixed-size array of 6 elements into a new dynamically allocated array.
+function copy6(int256[6] memory input) pure returns (int256[] memory) {
+    uint256 len = input.length;
+    int256[] memory output = new int256[](len);
+    for (uint256 i; i < len; i++) {
+        output[i] = input[i];
+    }
+    return output;
+}
+
+/// @dev Copies a fixed-size array of 7 elements into a new dynamically allocated array.
+function copy7(int256[7] memory input) pure returns (int256[] memory) {
+    uint256 len = input.length;
+    int256[] memory output = new int256[](len);
+    for (uint256 i; i < len; i++) {
+        output[i] = input[i];
+    }
+    return output;
+}
+
+/// @dev Copies a fixed-size array of 8 elements into a new dynamically allocated array.
+function copy8(int256[8] memory input) pure returns (int256[] memory) {
+    uint256 len = input.length;
+    int256[] memory output = new int256[](len);
+    for (uint256 i; i < len; i++) {
+        output[i] = input[i];
+    }
+    return output;
+}
+
+/// @dev Copies a fixed-size array of 9 elements into a new dynamically allocated array.
+function copy9(int256[9] memory input) pure returns (int256[] memory) {
+    uint256 len = input.length;
+    int256[] memory output = new int256[](len);
+    for (uint256 i; i < len; i++) {
+        output[i] = input[i];
+    }
+    return output;
+}
diff --git a/src/libraries/PolynomialLib.sol b/src/libraries/PolynomialLib.sol
new file mode 100644
index 0000000..091353e
--- /dev/null
+++ b/src/libraries/PolynomialLib.sol
@@ -0,0 +1,49 @@
+// SPDX-License-Identifier: AGPL-3.0-only
+pragma solidity ^0.8.4;
+
+import "./RayMathLib.sol";
+
+////////////////
+//  Errors    //
+////////////////
+
+error InvalidLength();
+
+////////////////
+// Functions  //
+////////////////
+
+/// @dev Evaluates a polynomial of the form:
+///      y = coef[0] + coef[1]*x + coef[2]*x^2 + ... + coef[N]*x^N
+///      The function is evaluated by Horner's method.
+/// @custom:reference Inspired by https://github.com/jeremybarnes/cephes/blob/master/cprob/polevl.c
+/// @param x The value at which the polynomial is evaluated.
+/// @param coef The coefficients of the polynomial.
+/// @param N The degree of the polynomial.
+/// @return y The value of the polynomial at x.
+function polevl(int256 x, int256[] memory coef, uint256 N) pure returns (int256 y) {
+    if (N >= coef.length) revert InvalidLength();
+
+    y = coef[0];
+    for (uint256 i = 1; i <= N; i++) {
+        y = mulfp(y, x) + coef[i];
+    }
+
+    return y;
+}
+
+/// @dev Evaluate polynomial when coefficient of x is 1.0. Otherwise same as polevl.
+/// @custom:reference Inspired by https://github.com/jeremybarnes/cephes/blob/master/cprob/polevl.c
+/// @param x The value at which the polynomial is evaluated.
+/// @param coef The coefficients of the polynomial.
+/// @param N The degree of the polynomial.
+/// @return y The value of the polynomial at x.
+function p1evl(int256 x, int256[] memory coef, uint256 N) pure returns (int256 y) {
+    if (N > coef.length) revert InvalidLength();
+
+    y = x + coef[0];
+    for (uint256 i = 1; i < N; i++) {
+        y = mulfp(y, x) + coef[i];
+    }
+    return y;
+}
diff --git a/src/libraries/RayMathLib.sol b/src/libraries/RayMathLib.sol
new file mode 100644
index 0000000..4156be4
--- /dev/null
+++ b/src/libraries/RayMathLib.sol
@@ -0,0 +1,64 @@
+// SPDX-License-Identifier: AGPL-3.0-only
+pragma solidity ^0.8.4;
+
+import "solmate/utils/FixedPointMathLib.sol";
+
+////////////////
+// Constants  //
+////////////////
+
+int256 constant RAY = 1e27;
+int256 constant RAY_HALF = 0.5e27;
+int256 constant RAY_ONE = RAY;
+int256 constant RAY_TWO = 2e27;
+int256 constant RAY_EIGHT = 8e27;
+int256 constant RAY_SQRT2 = 1.41421356237309504880168872e27; // sqrt(2)
+int256 constant RAY_SQRT2PI = 2.50662827463100050242e27; // sqrt(2pi)
+int256 constant RAY_MAXLOG = 7.09782712893383996843e29; // ln(2^255)
+
+////////////////
+// Arithmetic //
+////////////////
+
+/// @dev Fixed point multiplication, truncates.
+function mulfp(int256 a, int256 b) pure returns (int256) {
+    return (a * b) / RAY;
+}
+
+/// @dev Fixed point multiplication, truncates.
+function mulfp(uint256 a, uint256 b) pure returns (uint256) {
+    return (a * b) / uint256(RAY);
+}
+
+/// @dev Fixed point division, truncates.
+function divfp(int256 a, int256 b) pure returns (int256) {
+    return a * RAY / b;
+}
+
+/// @dev Fixed point division, truncates.
+function divfp(uint256 a, uint256 b) pure returns (uint256) {
+    return a * uint256(RAY) / b;
+}
+
+////////////////
+// Functions  //
+////////////////
+
+/// @dev Returns the absolute value of a number.
+/// todo: do we need type checking?
+function absolute(int256 x) view returns (int256) {
+    if (x < 0) return -x;
+    else return x;
+}
+
+/// @dev Up to 1E18 precision.
+function sqrtfp(int256 x) view returns (int256) {
+    x = x / 1e9;
+    return int256(FixedPointMathLib.sqrt(uint256(x))) * 1e18; // todo: fix, this is temp
+}
+
+/// @dev todo
+function logfp(int256 x) view returns (int256) {
+    x = x / 1e9 + 1;
+    return int256(FixedPointMathLib.lnWad(x)) * 1e9; // todo: fix, this is temp
+}
diff --git a/src/test/Ndtr.t.sol b/src/test/Ndtr.t.sol
new file mode 100644
index 0000000..7cbca66
--- /dev/null
+++ b/src/test/Ndtr.t.sol
@@ -0,0 +1,220 @@
+// SPDX-License-Identifier: GPL-3.0-only
+pragma solidity ^0.8.4;
+
+import "forge-std/Test.sol";
+import "../Ndtr.sol";
+
+import "../Gaussian.sol";
+
+/// @dev Hardcoded expected values computed manually using `normalcdlower` operation in https://keisan.casio.com/calculator
+contract TestNdtr is Test {
+    int256 constant NDTR_ACCURACY = 1e18;
+
+    /// @notice Compares two in256 values up to a precision with a base of RAY.
+    function assertEqPrecision(int256 a, int256 b, int256 precision, string memory message) internal {
+        // Gets the digits passed the precision end point.
+        uint256 remainder0 = mulmod(uint256(a), uint256(precision), uint256(RAY));
+        uint256 remainder1 = mulmod(uint256(b), uint256(precision), uint256(RAY));
+
+        // Add one to the remainder to round up, in the case it is 99...99.
+        remainder0++;
+        remainder1++;
+
+        // Converts units to precision.
+        a = a * precision / RAY;
+        b = b * precision / RAY;
+
+        // Rounds up if remainder is >= 0.5.
+        if (int256(remainder0) >= RAY_HALF) a++;
+        if (int256(remainder1) >= RAY_HALF) b++;
+
+        assertEq(a, b, message);
+    }
+
+    function test_ndtr() public {
+        assertEqPrecision(Ndtr.ndtr(RAY), int256(0.841344746068542948585232545e27), NDTR_ACCURACY, "ndtr-1");
+        assertEqPrecision(Ndtr.ndtr(RAY * 15 / 10), int256(0.933192798731141934488289354e27), NDTR_ACCURACY, "ndtr-1.5");
+
+        assertEqPrecision(Ndtr.ndtr(RAY * 2), int256(0.977249868051820792947945109e27), NDTR_ACCURACY, "ndtr-2");
+        assertEqPrecision(Ndtr.ndtr(RAY * 10), int256(0.999999999999999999999999999e27), NDTR_ACCURACY, "ndtr-10");
+        assertEqPrecision(
+            Ndtr.ndtr(RAY / 10000000), int256(0.500000039894228040143796598e27), NDTR_ACCURACY, "ndtr-0.0000001"
+        );
+        assertEqPrecision(Ndtr.ndtr(RAY / 100), int256(0.503989356314631603924543868e27), NDTR_ACCURACY, "ndtr-0.01");
+        // negative
+        assertEqPrecision(
+            Ndtr.ndtr(-RAY / 10000000), int256(0.499999960105771959856203402e27), NDTR_ACCURACY, "ndtr--0.0000001"
+        );
+        assertEqPrecision(Ndtr.ndtr(-RAY / 100), int256(0.496010643685368396075456132e27), NDTR_ACCURACY, "ndtr--0.01");
+        // 0
+        assertEqPrecision(Ndtr.ndtr(0), int256(0.5e27), NDTR_ACCURACY, "ndtr-0");
+        // random values at scale of 1e27
+        assertEqPrecision(
+            Ndtr.ndtr(0.123456789e27), int256(0.5491273050781420888711e27), NDTR_ACCURACY, "ndtr-0.123456789"
+        );
+        assertEqPrecision(
+            Ndtr.ndtr(0.987654321e27), int256(0.83833901356624443490786e27), NDTR_ACCURACY, "ndtr-0.987654321"
+        );
+
+        assertEqPrecision(
+            Ndtr.ndtr(-0.123456789e27), int256(0.4508726949218579111288e27), NDTR_ACCURACY, "ndtr--0.123456789"
+        );
+        assertEqPrecision(
+            Ndtr.ndtr(-0.987654321e27), int256(0.16166098643375556509213e27), NDTR_ACCURACY, "ndtr--0.987654321"
+        );
+    }
+
+    /// @dev Tests the edge case where the input is an integer near zero on the ray scale
+    function test_ndtr_near_zero() public {
+        assertEqPrecision(Ndtr.ndtr(RAY / NDTR_ACCURACY), int256(0.5e27), NDTR_ACCURACY, "ndtr-1e-18");
+        assertEqPrecision(
+            Ndtr.ndtr(-RAY / NDTR_ACCURACY), int256(0.499999999999999999999999999e27), NDTR_ACCURACY, "ndtr--1e-18"
+        );
+
+        assertEqPrecision(Ndtr.ndtr(RAY / NDTR_ACCURACY + 1), int256(0.5e27), NDTR_ACCURACY, "ndtr-1e-18+1");
+        assertEqPrecision(
+            Ndtr.ndtr(-RAY / NDTR_ACCURACY - 1),
+            int256(0.499999999999999999999999999e27),
+            NDTR_ACCURACY,
+            "ndtr--1e-18-1"
+        );
+
+        assertEqPrecision(
+            Ndtr.ndtr(RAY / NDTR_ACCURACY * 1e1), int256(0.500000000000000003989422804e27), NDTR_ACCURACY, "ndtr-1e-17"
+        );
+        assertEqPrecision(
+            Ndtr.ndtr(-RAY / NDTR_ACCURACY * 1e1),
+            int256(0.499999999999999996010577195e27),
+            NDTR_ACCURACY,
+            "ndtr--1e-17"
+        );
+
+        assertEqPrecision(
+            Ndtr.ndtr(RAY / NDTR_ACCURACY * 1e2), int256(0.50000000000000003989422804e27), NDTR_ACCURACY, "ndtr-1e-16"
+        );
+        assertEqPrecision(
+            Ndtr.ndtr(-RAY / NDTR_ACCURACY * 1e2),
+            int256(0.499999999999999960105771959e27),
+            NDTR_ACCURACY,
+            "ndtr--1e-16"
+        );
+
+        assertEqPrecision(
+            Ndtr.ndtr(RAY / NDTR_ACCURACY * 1e3), int256(0.500000000000000398942280401e27), NDTR_ACCURACY, "ndtr-1e-15"
+        );
+        assertEqPrecision(
+            Ndtr.ndtr(-RAY / NDTR_ACCURACY * 1e3),
+            int256(0.499999999999999601057719598e27),
+            NDTR_ACCURACY,
+            "ndtr--1e-15"
+        );
+
+        assertEqPrecision(
+            Ndtr.ndtr(RAY / NDTR_ACCURACY * 1e4), int256(0.500000000000003989422804014e27), NDTR_ACCURACY, "ndtr-1e-14"
+        );
+        assertEqPrecision(
+            Ndtr.ndtr(-RAY / NDTR_ACCURACY * 1e4),
+            int256(0.499999999999996010577195985e27),
+            NDTR_ACCURACY,
+            "ndtr--1e-14"
+        );
+
+        assertEqPrecision(
+            Ndtr.ndtr(RAY / NDTR_ACCURACY * 1e5), int256(0.500000000000039894228040143e27), NDTR_ACCURACY, "ndtr-1e-13"
+        );
+
+        assertEqPrecision(
+            Ndtr.ndtr(-RAY / NDTR_ACCURACY * 1e5),
+            int256(0.499999999999960105771959856e27),
+            NDTR_ACCURACY,
+            "ndtr--1e-13"
+        );
+    }
+
+    function test_ndtr_far_zero() public {
+        assertEqPrecision(Ndtr.ndtr(RAY * 2), int256(0.977249868051820792947945109e27), NDTR_ACCURACY, "ndtr-2");
+        assertEqPrecision(Ndtr.ndtr(RAY * 3), int256(0.998650101968369905473348185e27), NDTR_ACCURACY, "ndtr-3");
+        assertEqPrecision(Ndtr.ndtr(RAY * 4), int256(0.999968328758166880078746229e27), NDTR_ACCURACY, "ndtr-4");
+        assertEqPrecision(Ndtr.ndtr(RAY * 5), int256(0.999999713348428120806088326e27), NDTR_ACCURACY, "ndtr-5");
+        assertEqPrecision(Ndtr.ndtr(RAY * 6), int256(0.999999999013412354962301859e27), NDTR_ACCURACY, "ndtr-6");
+        assertEqPrecision(Ndtr.ndtr(RAY * 7), int256(0.999999999998720187456114164e27), NDTR_ACCURACY, "ndtr-7");
+        assertEqPrecision(Ndtr.ndtr(RAY * 8), int256(0.999999999999999377903942572e27), NDTR_ACCURACY, "ndtr-8");
+        assertEqPrecision(Ndtr.ndtr(RAY * 9), int256(0.999999999999999999887141159e27), NDTR_ACCURACY, "ndtr-9");
+        assertEqPrecision(Ndtr.ndtr(RAY * 10), int256(0.999999999999999999887141159e27), NDTR_ACCURACY, "ndtr-10");
+
+        // Negative
+        assertEqPrecision(Ndtr.ndtr(-RAY * 2), int256(0.022750131948179207200282637e27), NDTR_ACCURACY, "ndtr--2");
+        assertEqPrecision(Ndtr.ndtr(-RAY * 3), int256(0.001349898031630094526651814e27), NDTR_ACCURACY, "ndtr--3");
+        assertEqPrecision(Ndtr.ndtr(-RAY * 4), int256(3.167124183311992125377e22), NDTR_ACCURACY, "ndtr--4");
+        assertEqPrecision(Ndtr.ndtr(-RAY * 5), int256(2.86651571879193911673e20), NDTR_ACCURACY, "ndtr--5");
+        assertEqPrecision(Ndtr.ndtr(-RAY * 6), int256(9.8658764503769814e17), NDTR_ACCURACY, "ndtr--6");
+        assertEqPrecision(Ndtr.ndtr(-RAY * 7), int256(1.279812543885835e15), NDTR_ACCURACY, "ndtr--7");
+        assertEqPrecision(Ndtr.ndtr(-RAY * 8), int256(6.22096057427e11), NDTR_ACCURACY, "ndtr--8");
+        assertEqPrecision(Ndtr.ndtr(-RAY * 9), int256(1.1285884e8), NDTR_ACCURACY, "ndtr--9");
+        assertEqPrecision(Ndtr.ndtr(-RAY * 10), int256(7.619e3), NDTR_ACCURACY, "ndtr--10");
+    }
+
+    /// @dev Expected `b` values computed manually using `normalcdlower` operation in https://keisan.casio.com/calculator
+    function test_ndtr_near_half() public {
+        assertEqPrecision(Ndtr.ndtr(RAY / 2), int256(0.69146246127401310363770461e27), NDTR_ACCURACY, "ndtr-0.5");
+        assertEqPrecision(
+            Ndtr.ndtr(RAY * 7 / 12), int256(0.720165536400294250547948904e27), NDTR_ACCURACY, "ndtr-0.583333"
+        );
+        assertEqPrecision(
+            Ndtr.ndtr(RAY * 2 / 3), int256(0.747507462453077086935938175e27), NDTR_ACCURACY, "ndtr-0.666666"
+        );
+        assertEqPrecision(
+            Ndtr.ndtr(RAY * 5 / 6), int256(0.7976716190363569746314664e27), NDTR_ACCURACY, "ndtr-0.833333"
+        );
+
+        assertEqPrecision(Ndtr.ndtr(-RAY / 2), int256(0.308537538725986896362295389e27), NDTR_ACCURACY, "ndtr--0.5");
+        assertEqPrecision(
+            Ndtr.ndtr(-RAY * 7 / 12), int256(0.279834463599705749452051095e27), NDTR_ACCURACY, "ndtr--0.583333"
+        );
+        assertEqPrecision(
+            Ndtr.ndtr(-RAY * 2 / 3), int256(0.252492537546922913064061824e27), NDTR_ACCURACY, "ndtr--0.666666"
+        );
+        assertEqPrecision(
+            Ndtr.ndtr(-RAY * 5 / 6), int256(0.202328380963643025368533599e27), NDTR_ACCURACY, "ndtr--0.833333"
+        );
+    }
+
+    int256 constant UPPER_NDTR_BOUND = 10 * RAY;
+    int256 constant LOWER_NDTR_BOUND = -10 * RAY;
+
+    /// note: fails strictly at values above abs(1), i.e. a == b
+    function testFuzz_ndtr_monotonically_increasing(int256 x) public {
+        x = bound(x, LOWER_NDTR_BOUND, UPPER_NDTR_BOUND);
+
+        int256 a = Ndtr.ndtr(x);
+        int256 b = Ndtr.ndtr(x + 1);
+        console.logInt(a);
+        console.logInt(b);
+
+        assertTrue(a <= b, "ndtr-monotonically-increasing"); // note: not strict?
+    }
+
+    /// note: fails strictly at values above abs(1), i.e. a == b
+    function testFuzz_ndtr_monotonically_decreasing(int256 x) public {
+        x = bound(x, LOWER_NDTR_BOUND, UPPER_NDTR_BOUND);
+
+        int256 a = Ndtr.ndtr(x);
+        int256 b = Ndtr.ndtr(x - 1);
+        console.logInt(a);
+        console.logInt(b);
+
+        assertTrue(a >= b, "ndtr-monotonically-decreasing"); // note: not strict?
+    }
+
+    /// todo: update this test... it will fail because ndtr is better than the reference, so not equal!
+    function testFuzz_ndtr_reference(int256 x) public {
+        x = bound(x, -int256(2 * 1e27), int256(2 * 1e27));
+        int256 y0 = Ndtr.ndtr(x);
+        x /= 1e9;
+        int256 y1 = Gaussian.cdf(x);
+        console.logInt(y0);
+        console.logInt(y1);
+
+        assertEqPrecision(y0, y1 * 1e9, 1e9, "ndtr-reference");
+    }
+}
diff --git a/src/test/Ndtri.t.sol b/src/test/Ndtri.t.sol
new file mode 100644
index 0000000..d2a7209
--- /dev/null
+++ b/src/test/Ndtri.t.sol
@@ -0,0 +1,145 @@
+// SPDX-License-Identifier: GPL-3.0-only
+pragma solidity ^0.8.4;
+
+import "forge-std/Test.sol";
+import "../Ndtri.sol";
+
+/// @dev Hardcoded expected values computed manually using `normalcdlower` operation in https://keisan.casio.com/calculator
+contract TestNdtri is Test {
+    int256 constant NDTRI_ACCURACY = 1e18;
+
+    /// @notice Compares two in256 values up to a precision with a base of RAY.
+    function assertEqPrecision(int256 a, int256 b, int256 precision, string memory message) internal {
+        assertEq(a * precision / RAY, b * precision / RAY, message);
+    }
+
+    /// @dev Tests various inputs on the domain [0, 1], where input is a probability between 0 and 1.
+    function test_ndtri() public {
+        assertEqPrecision(Ndtri.ndtri(RAY / 2), int256(0), NDTRI_ACCURACY, "ndtri-0.5");
+        assertEqPrecision(Ndtri.ndtri(RAY / 4), int256(-0.674489750196081743202227014e27), NDTRI_ACCURACY, "ndtri-0.25");
+        assertEqPrecision(Ndtri.ndtri(RAY / 8), int256(-1.15034938037600817829676531e27), NDTRI_ACCURACY, "ndtri-0.125");
+        assertEqPrecision(
+            Ndtri.ndtri(RAY / 16), int256(-1.53412054435254631170839905e27), NDTRI_ACCURACY, "ndtri-0.0625"
+        );
+        assertEqPrecision(
+            Ndtri.ndtri(RAY * 77 / 93), int256(0.946122671364591906727098902e27), NDTRI_ACCURACY, "ndtri-0.827956"
+        );
+        assertEqPrecision(
+            Ndtri.ndtri(RAY * 11 / 85), int256(-1.1291761577077979287868872e27), NDTRI_ACCURACY, "ndtri-0.129412"
+        );
+    }
+
+    /// Starting at 1e27, decrements by atleast RAY / NDTRI_ACCURACY = 1e9, since thats the lowest unit
+    /// that has precision.
+    function test_ndtri_near_one() public {
+        int256 decrement = RAY / NDTRI_ACCURACY;
+        assertEqPrecision(
+            Ndtri.ndtri(RAY - decrement),
+            int256(8.75729034878231506388112862e27),
+            NDTRI_ACCURACY,
+            "ndtri-0.999999999999999999000000000"
+        );
+        assertEqPrecision(
+            Ndtri.ndtri(RAY - decrement * 1e1), int256(8.49379322410959807444471881e27), NDTRI_ACCURACY, "ndtri-1-1e-17"
+        );
+        assertEqPrecision(
+            Ndtri.ndtri(RAY - decrement * 1e2), int256(8.22208221613043561267585878e27), NDTRI_ACCURACY, "ndtri-1-1e-16"
+        );
+        assertEqPrecision(
+            Ndtri.ndtri(RAY - decrement * 1e3), int256(7.94134532617099678096674357e27), NDTRI_ACCURACY, "ndtri-1-1e-15"
+        );
+        assertEqPrecision(
+            Ndtri.ndtri(RAY - decrement * 1e4), int256(7.65062809293526881641896581e27), NDTRI_ACCURACY, "ndtri-1-1e-14"
+        );
+        assertEqPrecision(
+            Ndtri.ndtri(RAY - decrement * 1e5), int256(7.34879610280067751753910726e27), NDTRI_ACCURACY, "ndtri-1-1e-13"
+        );
+    }
+
+    /// @dev Tests the edge case where the input is an integer near zero on the ray scale
+    function test_ndtri_near_zero() public {
+        assertEqPrecision(
+            Ndtri.ndtri(RAY / NDTRI_ACCURACY + 1), -int256(8.7572903487823e27), NDTRI_ACCURACY, "ndtri-1e-18"
+        );
+        assertEqPrecision(
+            Ndtri.ndtri(RAY / NDTRI_ACCURACY * 1e1), -int256(8.493793224109598e27), NDTRI_ACCURACY, "ndtri-1e-17"
+        );
+        assertEqPrecision(
+            Ndtri.ndtri(RAY / NDTRI_ACCURACY * 1e2), -int256(8.2220822161304356e27), NDTRI_ACCURACY, "ndtri-1e-16"
+        );
+        assertEqPrecision(
+            Ndtri.ndtri(RAY / NDTRI_ACCURACY * 1e3), -int256(7.9413453261709968e27), NDTRI_ACCURACY, "ndtri-1e-15"
+        );
+        assertEqPrecision(
+            Ndtri.ndtri(RAY / NDTRI_ACCURACY * 1e4), -int256(7.65062809293526882e27), NDTRI_ACCURACY, "ndtri-1e-14"
+        );
+        assertEqPrecision(
+            Ndtri.ndtri(RAY / NDTRI_ACCURACY * 1e5), -int256(7.348796102800677518e27), NDTRI_ACCURACY, "ndtri-1e-13"
+        );
+    }
+
+    function test_ndtri_far_zero() public {}
+
+    /// @dev Expected `b` values computed manually using `normalcdlower` operation in https://keisan.casio.com/calculator
+    function test_ndtri_near_half() public {
+        assertEqPrecision(
+            Ndtri.ndtri(RAY * 3 / 5), int256(0.253347103135799798798196181e27), NDTRI_ACCURACY, "ndtri-0.666666"
+        );
+        assertEqPrecision(
+            Ndtri.ndtri(RAY * 7 / 12), int256(0.210428394247924723324540643e27), NDTRI_ACCURACY, "ndtri-0.583333"
+        );
+        assertEqPrecision(
+            Ndtri.ndtri(RAY * 2 / 3), int256(0.43072729929545749020594039e27), NDTRI_ACCURACY, "ndtri-0.666666"
+        );
+        assertEqPrecision(
+            Ndtri.ndtri(RAY * 5 / 6), int256(0.96742156610170103955040122e27), NDTRI_ACCURACY, "ndtri-0.833333"
+        );
+    }
+
+    int256 constant UPPER_NDTRI_BOUND = RAY - 1 - 1 - 1; // Bound is RAY - 1, but we want to increase by 1 so need at least 2
+    int256 constant LOWER_NDTRI_BOUND = 1 + 1 + 1; // Bound is 1, but we want to decrease by 1, so need at least 2
+
+    /// note: fails strictly at values above abs(1), i.e. a == b
+    function testFuzz_ndtri_monotonically_increasing(int256 x) public {
+        x = bound(x, LOWER_NDTRI_BOUND, UPPER_NDTRI_BOUND);
+
+        int256 a = Ndtri.ndtri(x);
+        int256 b = Ndtri.ndtri(x + 1);
+        console.logInt(a);
+        console.logInt(b);
+
+        assertTrue(a <= b, "ndtri-monotonically-increasing"); // note: not strict?
+    }
+
+    /// note: fails strictly at values above abs(1), i.e. a == b
+    function testFuzz_ndtri_monotonically_decreasing(int256 x) public {
+        x = bound(x, LOWER_NDTRI_BOUND, UPPER_NDTRI_BOUND);
+
+        int256 a = Ndtri.ndtri(x);
+        int256 b = Ndtri.ndtri(x - 1);
+        console.logInt(a);
+        console.logInt(b);
+
+        assertTrue(a >= b, "ndtri-monotonically-decreasing"); // note: not strict?
+    }
+
+    function test_ndtri_equals_one_reverts() public {
+        vm.expectRevert(Ndtri.MaxNumError.selector);
+        Ndtri.ndtri(RAY);
+    }
+
+    function test_ndtri_gte_one_reverts() public {
+        vm.expectRevert(Ndtri.MaxNumError.selector);
+        Ndtri.ndtri(RAY + 1);
+    }
+
+    function test_ndtri_equals_zero_reverts() public {
+        vm.expectRevert(Ndtri.MaxNumError.selector);
+        Ndtri.ndtri(0);
+    }
+
+    function test_ndtri_lte_zero_reverts() public {
+        vm.expectRevert(Ndtri.MaxNumError.selector);
+        Ndtri.ndtri(0 - 1);
+    }
+}
diff --git a/src/test/RayMathLib.t.sol b/src/test/RayMathLib.t.sol
new file mode 100644
index 0000000..24c2b78
--- /dev/null
+++ b/src/test/RayMathLib.t.sol
@@ -0,0 +1,15 @@
+// SPDX-License-Identifier: GPL-3.0-only
+pragma solidity ^0.8.4;
+
+import "forge-std/Test.sol";
+import "solmate/utils/FixedPointMathLib.sol";
+import "../libraries/RayMathLib.sol" as RayMathLib;
+
+contract TestRayMathLib is Test {
+    function test_logfp() public {
+        int256 x = 3.5e27;
+        int256 y = RayMathLib.logfp(x);
+        console.logInt(FixedPointMathLib.lnWad(3.5e18));
+        assertEq(y, 1.25276296849536799568812062e27, "logfp-3.5e27");
+    }
+}

From cd55a0c0c6da290d43bc834753359b17f933fe97 Mon Sep 17 00:00:00 2001
From: alex <alexangel.j@gmail.com>
Date: Sun, 21 May 2023 16:08:22 -0700
Subject: [PATCH 2/9] feat(monotonicity-tests): adds invariant monotonicity
 tests

---
 foundry.toml                     |   2 +-
 src/libraries/RayMathLib.sol     |  96 +++++++++++++++++++++++++-
 src/test/Ndtr.t.sol              |  35 ++++++++++
 src/test/RMM01Monotonicity.t.sol | 113 +++++++++++++++++++++++++++++++
 src/test/RayMathLib.t.sol        |   9 +++
 5 files changed, 252 insertions(+), 3 deletions(-)
 create mode 100644 src/test/RMM01Monotonicity.t.sol

diff --git a/foundry.toml b/foundry.toml
index 06d6fdb..1d0bc4b 100644
--- a/foundry.toml
+++ b/foundry.toml
@@ -11,5 +11,5 @@ remappings = [
 via_ir = false
 
 [fuzz]
-runs = 100000
+runs = 256
 max_test_rejects = 4000000
diff --git a/src/libraries/RayMathLib.sol b/src/libraries/RayMathLib.sol
index 4156be4..dc7c466 100644
--- a/src/libraries/RayMathLib.sol
+++ b/src/libraries/RayMathLib.sol
@@ -53,8 +53,10 @@ function absolute(int256 x) view returns (int256) {
 
 /// @dev Up to 1E18 precision.
 function sqrtfp(int256 x) view returns (int256) {
-    x = x / 1e9;
-    return int256(FixedPointMathLib.sqrt(uint256(x))) * 1e18; // todo: fix, this is temp
+    uint256 result = sqrt(uint256(x));
+    if (result > uint256(type(int256).max)) revert("overflow");
+
+    return int256(result) * 1e14; // Multiplies by 1e14 since the sqrt of 1e27 units is 1e13.5.
 }
 
 /// @dev todo
@@ -62,3 +64,93 @@ function logfp(int256 x) view returns (int256) {
     x = x / 1e9 + 1;
     return int256(FixedPointMathLib.lnWad(x)) * 1e9; // todo: fix, this is temp
 }
+
+/// @notice Calculates the square root of x using the Babylonian method.
+///
+/// @dev See https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Babylonian_method.
+///
+/// Notes:
+/// - If x is not a perfect square, the result is rounded down.
+/// - Credits to OpenZeppelin for the explanations in comments below.
+///
+/// @param x The uint256 number for which to calculate the square root.
+/// @return result The result as a uint256.
+/// @custom:smtchecker abstract-function-nondet
+function sqrt(uint256 x) pure returns (uint256 result) {
+    if (x == 0) {
+        return 0;
+    }
+
+    // For our first guess, we calculate the biggest power of 2 which is smaller than the square root of x.
+    //
+    // We know that the "msb" (most significant bit) of x is a power of 2 such that we have:
+    //
+    // $$
+    // msb(x) <= x <= 2*msb(x)$
+    // $$
+    //
+    // We write $msb(x)$ as $2^k$, and we get:
+    //
+    // $$
+    // k = log_2(x)
+    // $$
+    //
+    // Thus, we can write the initial inequality as:
+    //
+    // $$
+    // 2^{log_2(x)} <= x <= 2*2^{log_2(x)+1} \\
+    // sqrt(2^k) <= sqrt(x) < sqrt(2^{k+1}) \\
+    // 2^{k/2} <= sqrt(x) < 2^{(k+1)/2} <= 2^{(k/2)+1}
+    // $$
+    //
+    // Consequently, $2^{log_2(x) /2} is a good first approximation of sqrt(x) with at least one correct bit.
+    uint256 xAux = uint256(x);
+    result = 1;
+    if (xAux >= 2 ** 128) {
+        xAux >>= 128;
+        result <<= 64;
+    }
+    if (xAux >= 2 ** 64) {
+        xAux >>= 64;
+        result <<= 32;
+    }
+    if (xAux >= 2 ** 32) {
+        xAux >>= 32;
+        result <<= 16;
+    }
+    if (xAux >= 2 ** 16) {
+        xAux >>= 16;
+        result <<= 8;
+    }
+    if (xAux >= 2 ** 8) {
+        xAux >>= 8;
+        result <<= 4;
+    }
+    if (xAux >= 2 ** 4) {
+        xAux >>= 4;
+        result <<= 2;
+    }
+    if (xAux >= 2 ** 2) {
+        result <<= 1;
+    }
+
+    // At this point, `result` is an estimation with at least one bit of precision. We know the true value has at
+    // most 128 bits, since it is the square root of a uint256. Newton's method converges quadratically (precision
+    // doubles at every iteration). We thus need at most 7 iteration to turn our partial result with one bit of
+    // precision into the expected uint128 result.
+    unchecked {
+        result = (result + x / result) >> 1;
+        result = (result + x / result) >> 1;
+        result = (result + x / result) >> 1;
+        result = (result + x / result) >> 1;
+        result = (result + x / result) >> 1;
+        result = (result + x / result) >> 1;
+        result = (result + x / result) >> 1;
+
+        // If x is not a perfect square, round the result toward zero.
+        uint256 roundedResult = x / result;
+        if (result >= roundedResult) {
+            result = roundedResult;
+        }
+    }
+}
diff --git a/src/test/Ndtr.t.sol b/src/test/Ndtr.t.sol
index 7cbca66..41b15fa 100644
--- a/src/test/Ndtr.t.sol
+++ b/src/test/Ndtr.t.sol
@@ -217,4 +217,39 @@ contract TestNdtr is Test {
 
         assertEqPrecision(y0, y1 * 1e9, 1e9, "ndtr-reference");
     }
+
+    function testFuzz_erfc_ReturnsTwoWhenInputIsTooLow(int256 x) public {
+        vm.assume(x <= Gaussian.ERFC_DOMAIN_LOWER);
+
+        // TODO: Investigate why the error selector is 0x4d2d75b1
+        // instead of 0x35278d12
+        int256 y = Gaussian.erfc(x);
+        assertEq(y, int256(2 ether), "erfc-not-two");
+    }
+
+    function testFuzz_erfc_ReturnsZeroWhenInputIsTooHigh(int256 x) public {
+        vm.assume(x >= Gaussian.ERFC_DOMAIN_UPPER);
+        int256 y = Gaussian.erfc(x);
+        assertEq(y, 0, "erfc-not-zero");
+    }
+
+    function testFuzz_erfc_NegativeInputIsBounded(int256 x) public {
+        vm.assume(x > -1999999999999999998000000000000000002);
+        vm.assume(x < -0.0000001 ether);
+        int256 y = Gaussian.erfc(x);
+        assertGe(y, 1 ether);
+        assertLe(y, 2 ether);
+    }
+
+    function testFuzz_erfc_PositiveInputIsBounded(int256 x) public {
+        vm.assume(x > 0.0000001 ether);
+        vm.assume(x < 1999999999999999998000000000000000002);
+        int256 y = Gaussian.erfc(x);
+        assertGe(y, 0 ether);
+        assertLe(y, 1 ether);
+    }
+
+    function test_erfc_zero_input_returns_one() public {
+        assertEq(Gaussian.erfc(0), 1 ether);
+    }
 }
diff --git a/src/test/RMM01Monotonicity.t.sol b/src/test/RMM01Monotonicity.t.sol
new file mode 100644
index 0000000..a2aae65
--- /dev/null
+++ b/src/test/RMM01Monotonicity.t.sol
@@ -0,0 +1,113 @@
+// SPDX-License-Identifier: GPL-3.0-only
+pragma solidity ^0.8.4;
+
+import "forge-std/Test.sol";
+import "../Ndtr.sol";
+import "../Ndtri.sol";
+
+library RMM01 {
+    uint256 constant PRECISION = 1e27;
+    int256 internal constant YEAR = 31556952;
+
+    /// @dev Input x reserve is out of the bounds [0, 1].
+    error InvalidX();
+
+    /// @dev Up to 1E18 precision.
+    function getY(
+        uint256 x,
+        uint256 strikePriceRay,
+        uint256 volatilityRay,
+        uint256 timeRemainingSeconds,
+        int256 currentInvariant
+    ) internal view returns (uint256 y) {
+        if (x > PRECISION) revert InvalidX();
+
+        if (x == PRECISION) return uint256(currentInvariant);
+        if (x == 0) return uint256(int256(strikePriceRay) + currentInvariant);
+
+        if (timeRemainingSeconds == 0) {
+            uint256 product = mulfp(strikePriceRay, PRECISION - x);
+            y = uint256(int256(product) + currentInvariant);
+        } else {
+            int256 tauSeconds = divfp(int256(timeRemainingSeconds), YEAR);
+            int256 sqrtTau = sqrtfp(tauSeconds); // * 1e18; // sqrt(1e27) gives us a valuew with 1e9 units... so mulfp by 1e18 to get 1e27 units
+            logger.logInt(sqrtTau);
+            int256 sqrtTauVol = mulfp(int256(volatilityRay), sqrtTau);
+            logger.logInt(sqrtTauVol);
+            int256 phi = int256(PRECISION) - int256(x);
+            logger.logInt(phi);
+            int256 inverseCdf = Ndtri.ndtri(phi);
+            logger.logInt(inverseCdf);
+            int256 input = inverseCdf - sqrtTauVol;
+            logger.logInt(input);
+            int256 cdf = Ndtr.ndtr(input);
+            logger.logInt(cdf);
+            int256 product = mulfp(int256(strikePriceRay), cdf);
+            y = uint256(product + currentInvariant);
+        }
+
+        return y;
+    }
+
+    /// @dev Computes `k` in `k = y - KΦ(Φ⁻¹(1-x) - σ√τ)`.
+    function invariant(
+        uint256 x,
+        uint256 y,
+        uint256 strikePriceRay,
+        uint256 volatilityRay,
+        uint256 timeRemainingSeconds,
+        int256 currentInvariant
+    ) internal view returns (int256) {
+        return int256(y) - int256(getY(x, strikePriceRay, volatilityRay, timeRemainingSeconds, currentInvariant));
+    }
+}
+
+contract TestRMM01Monotonicity is Test {
+    /// As x -> 1, y -> 0
+    function testFuzz_rmm_montonically_increasing(uint256 x, uint256 v, uint256 t) public {
+        x = bound(x, 2, 1e27 - 2); // ray
+        v = bound(v, 1, 1e7); // bps
+        t = bound(t, 1, 500 days); // seconds
+
+        // scales v up to ray since we bound by a small range in units of bps, since this is a realistic range
+        v = v * 1e27 / 1e4; // Divides by BPS units which is 1e4, since 1e4 == 100% == 1.0 == 1e27 ray.
+
+        uint256 strikePriceRay = 1e27; // simple strike price since the domain is large.
+        int256 k = 0;
+
+        uint256 y = RMM01.getY(x, strikePriceRay, v, t, k);
+        uint256 y2 = RMM01.getY(x + 1, strikePriceRay, v, t, k);
+        vm.assume(y > 0);
+
+        console.log("x", x);
+        console.log("y", y);
+        console.log("x2", x + 1);
+        console.log("y2", y2);
+
+        assertTrue(y <= y2, "y <= y2");
+    }
+
+    /// As x -> 0, y -> K
+    function testFuzz_rmm_montonically_decreasing(uint256 x, uint256 v, uint256 t) public {
+        x = bound(x, 2, 1e27 - 2); // ray
+        v = bound(v, 1, 1e7); // bps
+        t = bound(t, 1, 500 days); // seconds
+
+        // scales v up to ray since we bound by a small range in units of bps, since this is a realistic range
+        v = v * 1e27 / 1e4; // Divides by BPS units which is 1e4, since 1e4 == 100% == 1.0 == 1e27 ray.
+
+        uint256 strikePriceRay = 1e27; // simple strike price since the domain is large.
+        int256 k = 0;
+
+        uint256 y = RMM01.getY(x, strikePriceRay, v, t, k);
+        uint256 y2 = RMM01.getY(x - 1, strikePriceRay, v, t, k);
+        vm.assume(y > 0);
+
+        console.log("x", x);
+        console.log("y", y);
+        console.log("x2", x - 1);
+        console.log("y2", y2);
+
+        assertTrue(y2 <= y, "y2 <= y");
+    }
+}
diff --git a/src/test/RayMathLib.t.sol b/src/test/RayMathLib.t.sol
index 24c2b78..0e78b4c 100644
--- a/src/test/RayMathLib.t.sol
+++ b/src/test/RayMathLib.t.sol
@@ -12,4 +12,13 @@ contract TestRayMathLib is Test {
         console.logInt(FixedPointMathLib.lnWad(3.5e18));
         assertEq(y, 1.25276296849536799568812062e27, "logfp-3.5e27");
     }
+
+    function testFuzz_sqrtfp(int256 x) public {
+        vm.assume(x > 0);
+
+        int256 y = RayMathLib.sqrtfp(x); // Units of 1e14
+        int256 y2 = int256(FixedPointMathLib.sqrt(uint256(x))); // Units of 1e9
+
+        assertEq(y / 1e14, y2, "sqrtfp");
+    }
 }

From 2d697a5130b8eef23bfe487793a7248c7306c417 Mon Sep 17 00:00:00 2001
From: alex <alexangel.j@gmail.com>
Date: Sun, 21 May 2023 21:52:46 -0700
Subject: [PATCH 3/9] fix(sqrt): uses fixedpointmathlib sqrt, even though it
 has 1e9 precision

---
 cli/ndtr.ts                  | 35 +++++++++++++
 src/Ndtri.sol                |  7 ++-
 src/libraries/RayMathLib.sol | 99 +++---------------------------------
 src/test/Ndtri.t.sol         | 64 +++++++++++++++++++++--
 4 files changed, 106 insertions(+), 99 deletions(-)
 create mode 100644 cli/ndtr.ts

diff --git a/cli/ndtr.ts b/cli/ndtr.ts
new file mode 100644
index 0000000..16cd3f1
--- /dev/null
+++ b/cli/ndtr.ts
@@ -0,0 +1,35 @@
+async function main() {
+  const args = process.argv
+  const [, , flag, value] = args
+
+  if (value) {
+    const y0 = parseFloat(value)
+    ndtr(y0)
+  } else {
+    console.log('No value supplied...')
+  }
+}
+
+function ndtr(y0: number) {
+  let x, x0, z
+  let code = 1
+  let y = y0
+
+  if (y > 1.0 - 0.13533528323661269189) {
+    /* 0.135... = exp(-2) */
+    y = 1.0 - y
+    code = 0
+  }
+
+  let logy = Math.log(y)
+  x = Math.sqrt(-2.0 * Math.log(y))
+  x0 = x - Math.log(x) / x
+  z = 1.0 / x
+
+  console.log({ y0, y, logy, x, x0, z })
+}
+
+main().catch((err) => {
+  console.error(err)
+  process.exit(1)
+})
diff --git a/src/Ndtri.sol b/src/Ndtri.sol
index ceb73e3..db213ba 100644
--- a/src/Ndtri.sol
+++ b/src/Ndtri.sol
@@ -138,6 +138,8 @@ library Ndtri {
         int256 code = 1;
 
         y = y0;
+
+        // Use a different approximation if y > exp(-2).
         if (y > RAY_ONE - 0.13533528323661269189e27) {
             // 0.135... = exp(-2)
             y = RAY_ONE - y;
@@ -160,6 +162,9 @@ library Ndtri {
         // It also helps ensure that x will be a positive number because we're taking the square root.
         // The -2.0 * log(y) essentially gives the quantile of the exponential distribution,
         // which is then square rooted to approximate the quantile for the normal distribution.
+
+        // The precision of the sqrt output affects every value after this point.
+        // i.e. x0 will have precision up to the sqrt precision, and so on.
         x = sqrtfp(mulfp(-RAY_TWO, logfp(y))); // √(-2 * ln(y))
 
         // Part of a numerical approximation method (like a variant of Newton's method)
@@ -177,7 +182,7 @@ library Ndtri {
             x1 = z * polevl(z, P2_ARRAY, uint256(8)) / p1evl(z, Q2_ARRAY, uint256(8));
         }
 
-        x = x0 - y;
+        x = x0 - x1;
         if (code != 0) {
             x = -x;
         }
diff --git a/src/libraries/RayMathLib.sol b/src/libraries/RayMathLib.sol
index dc7c466..fa8b8f8 100644
--- a/src/libraries/RayMathLib.sol
+++ b/src/libraries/RayMathLib.sol
@@ -52,11 +52,14 @@ function absolute(int256 x) view returns (int256) {
 }
 
 /// @dev Up to 1E18 precision.
+/// todo: need more precise sqrt for 1e27 precision.
 function sqrtfp(int256 x) view returns (int256) {
-    uint256 result = sqrt(uint256(x));
-    if (result > uint256(type(int256).max)) revert("overflow");
+    x = x / 1e9; // Scale from 1e27 -> 1e18.
 
-    return int256(result) * 1e14; // Multiplies by 1e14 since the sqrt of 1e27 units is 1e13.5.
+    uint256 result = FixedPointMathLib.sqrt(uint256(x));
+    if (result > uint256(type(int256).max)) revert("sqrt-overflow");
+
+    return int256(result) * 1e18; // Multiplies by 1e18 since the sqrt of 1e18 units is 9, so multiply by 1e18 to get 1e27 units.
 }
 
 /// @dev todo
@@ -64,93 +67,3 @@ function logfp(int256 x) view returns (int256) {
     x = x / 1e9 + 1;
     return int256(FixedPointMathLib.lnWad(x)) * 1e9; // todo: fix, this is temp
 }
-
-/// @notice Calculates the square root of x using the Babylonian method.
-///
-/// @dev See https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Babylonian_method.
-///
-/// Notes:
-/// - If x is not a perfect square, the result is rounded down.
-/// - Credits to OpenZeppelin for the explanations in comments below.
-///
-/// @param x The uint256 number for which to calculate the square root.
-/// @return result The result as a uint256.
-/// @custom:smtchecker abstract-function-nondet
-function sqrt(uint256 x) pure returns (uint256 result) {
-    if (x == 0) {
-        return 0;
-    }
-
-    // For our first guess, we calculate the biggest power of 2 which is smaller than the square root of x.
-    //
-    // We know that the "msb" (most significant bit) of x is a power of 2 such that we have:
-    //
-    // $$
-    // msb(x) <= x <= 2*msb(x)$
-    // $$
-    //
-    // We write $msb(x)$ as $2^k$, and we get:
-    //
-    // $$
-    // k = log_2(x)
-    // $$
-    //
-    // Thus, we can write the initial inequality as:
-    //
-    // $$
-    // 2^{log_2(x)} <= x <= 2*2^{log_2(x)+1} \\
-    // sqrt(2^k) <= sqrt(x) < sqrt(2^{k+1}) \\
-    // 2^{k/2} <= sqrt(x) < 2^{(k+1)/2} <= 2^{(k/2)+1}
-    // $$
-    //
-    // Consequently, $2^{log_2(x) /2} is a good first approximation of sqrt(x) with at least one correct bit.
-    uint256 xAux = uint256(x);
-    result = 1;
-    if (xAux >= 2 ** 128) {
-        xAux >>= 128;
-        result <<= 64;
-    }
-    if (xAux >= 2 ** 64) {
-        xAux >>= 64;
-        result <<= 32;
-    }
-    if (xAux >= 2 ** 32) {
-        xAux >>= 32;
-        result <<= 16;
-    }
-    if (xAux >= 2 ** 16) {
-        xAux >>= 16;
-        result <<= 8;
-    }
-    if (xAux >= 2 ** 8) {
-        xAux >>= 8;
-        result <<= 4;
-    }
-    if (xAux >= 2 ** 4) {
-        xAux >>= 4;
-        result <<= 2;
-    }
-    if (xAux >= 2 ** 2) {
-        result <<= 1;
-    }
-
-    // At this point, `result` is an estimation with at least one bit of precision. We know the true value has at
-    // most 128 bits, since it is the square root of a uint256. Newton's method converges quadratically (precision
-    // doubles at every iteration). We thus need at most 7 iteration to turn our partial result with one bit of
-    // precision into the expected uint128 result.
-    unchecked {
-        result = (result + x / result) >> 1;
-        result = (result + x / result) >> 1;
-        result = (result + x / result) >> 1;
-        result = (result + x / result) >> 1;
-        result = (result + x / result) >> 1;
-        result = (result + x / result) >> 1;
-        result = (result + x / result) >> 1;
-
-        // If x is not a perfect square, round the result toward zero.
-        uint256 roundedResult = x / result;
-        if (result >= roundedResult) {
-            result = roundedResult;
-        }
-    }
-}
diff --git a/src/test/Ndtri.t.sol b/src/test/Ndtri.t.sol
index d2a7209..08e0f55 100644
--- a/src/test/Ndtri.t.sol
+++ b/src/test/Ndtri.t.sol
@@ -3,6 +3,7 @@ pragma solidity ^0.8.4;
 
 import "forge-std/Test.sol";
 import "../Ndtri.sol";
+import "../Gaussian.sol";
 
 /// @dev Hardcoded expected values computed manually using `normalcdlower` operation in https://keisan.casio.com/calculator
 contract TestNdtri is Test {
@@ -10,7 +11,23 @@ contract TestNdtri is Test {
 
     /// @notice Compares two in256 values up to a precision with a base of RAY.
     function assertEqPrecision(int256 a, int256 b, int256 precision, string memory message) internal {
-        assertEq(a * precision / RAY, b * precision / RAY, message);
+        // Gets the digits passed the precision end point.
+        uint256 remainder0 = mulmod(uint256(a), uint256(precision), uint256(RAY));
+        uint256 remainder1 = mulmod(uint256(b), uint256(precision), uint256(RAY));
+
+        // Add one to the remainder to round up, in the case it is 99...99.
+        remainder0++;
+        remainder1++;
+
+        // Converts units to precision.
+        a = a * precision / RAY;
+        b = b * precision / RAY;
+
+        // Rounds up if remainder is >= 0.5.
+        //if (int256(remainder0) >= RAY_HALF) a++;
+        //if (int256(remainder1) >= RAY_HALF) b++;
+
+        assertEq(a, b, message);
     }
 
     /// @dev Tests various inputs on the domain [0, 1], where input is a probability between 0 and 1.
@@ -27,12 +44,20 @@ contract TestNdtri is Test {
         assertEqPrecision(
             Ndtri.ndtri(RAY * 11 / 85), int256(-1.1291761577077979287868872e27), NDTRI_ACCURACY, "ndtri-0.129412"
         );
+
+        assertEqPrecision(
+            Ndtri.ndtri(729329433526774616220000000),
+            int256(0.610786196335546300516717471e27),
+            NDTRI_ACCURACY,
+            "ndtri-0.729329"
+        );
     }
 
     /// Starting at 1e27, decrements by atleast RAY / NDTRI_ACCURACY = 1e9, since thats the lowest unit
     /// that has precision.
     function test_ndtri_near_one() public {
         int256 decrement = RAY / NDTRI_ACCURACY;
+        console.logInt(RAY - decrement);
         assertEqPrecision(
             Ndtri.ndtri(RAY - decrement),
             int256(8.75729034878231506388112862e27),
@@ -99,7 +124,7 @@ contract TestNdtri is Test {
     int256 constant UPPER_NDTRI_BOUND = RAY - 1 - 1 - 1; // Bound is RAY - 1, but we want to increase by 1 so need at least 2
     int256 constant LOWER_NDTRI_BOUND = 1 + 1 + 1; // Bound is 1, but we want to decrease by 1, so need at least 2
 
-    /// note: fails strictly at values above abs(1), i.e. a == b
+    /// note: strict passes... does it pass with 10k fuzz runs though?
     function testFuzz_ndtri_monotonically_increasing(int256 x) public {
         x = bound(x, LOWER_NDTRI_BOUND, UPPER_NDTRI_BOUND);
 
@@ -108,10 +133,10 @@ contract TestNdtri is Test {
         console.logInt(a);
         console.logInt(b);
 
-        assertTrue(a <= b, "ndtri-monotonically-increasing"); // note: not strict?
+        assertTrue(a < b, "ndtri-monotonically-increasing");
     }
 
-    /// note: fails strictly at values above abs(1), i.e. a == b
+    /// note: strict passes... does it pass with 10k fuzz runs though?
     function testFuzz_ndtri_monotonically_decreasing(int256 x) public {
         x = bound(x, LOWER_NDTRI_BOUND, UPPER_NDTRI_BOUND);
 
@@ -120,7 +145,7 @@ contract TestNdtri is Test {
         console.logInt(a);
         console.logInt(b);
 
-        assertTrue(a >= b, "ndtri-monotonically-decreasing"); // note: not strict?
+        assertTrue(a > b, "ndtri-monotonically-decreasing");
     }
 
     function test_ndtri_equals_one_reverts() public {
@@ -142,4 +167,33 @@ contract TestNdtri is Test {
         vm.expectRevert(Ndtri.MaxNumError.selector);
         Ndtri.ndtri(0 - 1);
     }
+
+    function test_ndtri_lower() public {
+        assertEqPrecision(Ndtri.ndtri(0.1e27), int256(-1.28155156554460046696510332e27), NDTRI_ACCURACY, "ndtri-0.1");
+    }
+
+    /* /// todo: update this test... it will fail because ndtri is better than the reference, so not equal!
+    function testFuzz_ndtri_reference(int256 x) public {
+        x = bound(x, 2, 1e27 - 2);
+        int256 y0 = Ndtri.ndtri(x);
+        x /= 1e9;
+        int256 y1 = Gaussian.ppf(x);
+        console.logInt(y0);
+        console.logInt(y1);
+
+        assertEqPrecision(y0, y1 * 1e9, 1e9, "ndtri-reference");
+    }
+
+    /// todo: update this test... it will fail because ndtri is better than the reference, so not equal!
+    function testFuzz_ndtri_refrence_below_exp2(int256 x) public {
+        x = bound(x, 0.13533528323661269189e27 + 1, 1e27 - 0.13533528323661269189e27);
+
+        int256 y0 = Ndtri.ndtri(x);
+        x /= 1e9;
+        int256 y1 = Gaussian.ppf(x);
+        console.logInt(y0);
+        console.logInt(y1);
+
+        assertEqPrecision(y0, y1 * 1e9, 1e9, "ndtri-reference-below-exp2");
+    } */
 }

From 501d93d77f3e023c3de9c1a35ffae11d2dbf80d7 Mon Sep 17 00:00:00 2001
From: alex <alexangel.j@gmail.com>
Date: Sun, 21 May 2023 22:45:05 -0700
Subject: [PATCH 4/9] test(cleanup): cleans up some code and tests

---
 src/Ndtr.sol                     |   2 +-
 src/libraries/RayMathLib.sol     |  15 +++-
 src/test/Ndtr.t.sol              |  65 ++++++++--------
 src/test/RMM01Monotonicity.t.sol | 130 ++++++++++++++++++++++++-------
 4 files changed, 146 insertions(+), 66 deletions(-)

diff --git a/src/Ndtr.sol b/src/Ndtr.sol
index 67c076f..dc125b4 100644
--- a/src/Ndtr.sol
+++ b/src/Ndtr.sol
@@ -151,7 +151,7 @@ library Ndtr {
             else return 0;
         }
 
-        z = FixedPointMathLib.expWad(z / 1e9) * 1e9;
+        z = expfp(z);
 
         if (x < RAY_EIGHT) {
             int256[] memory input0 = copy9(abi.decode(P, (int256[9])));
diff --git a/src/libraries/RayMathLib.sol b/src/libraries/RayMathLib.sol
index fa8b8f8..1ff86b0 100644
--- a/src/libraries/RayMathLib.sol
+++ b/src/libraries/RayMathLib.sol
@@ -46,14 +46,14 @@ function divfp(uint256 a, uint256 b) pure returns (uint256) {
 
 /// @dev Returns the absolute value of a number.
 /// todo: do we need type checking?
-function absolute(int256 x) view returns (int256) {
+function absolute(int256 x) pure returns (int256) {
     if (x < 0) return -x;
     else return x;
 }
 
 /// @dev Up to 1E18 precision.
 /// todo: need more precise sqrt for 1e27 precision.
-function sqrtfp(int256 x) view returns (int256) {
+function sqrtfp(int256 x) pure returns (int256) {
     x = x / 1e9; // Scale from 1e27 -> 1e18.
 
     uint256 result = FixedPointMathLib.sqrt(uint256(x));
@@ -63,7 +63,16 @@ function sqrtfp(int256 x) view returns (int256) {
 }
 
 /// @dev todo
-function logfp(int256 x) view returns (int256) {
+function logfp(int256 x) pure returns (int256) {
     x = x / 1e9 + 1;
     return int256(FixedPointMathLib.lnWad(x)) * 1e9; // todo: fix, this is temp
 }
+
+/// @dev todo
+function expfp(int256 x) pure returns (int256) {
+    x = x / 1e9;
+
+    int256 result = FixedPointMathLib.expWad(x);
+
+    return result * 1e9; // todo: fix, this is temp, scale from 1E18 -> 1E27
+}
diff --git a/src/test/Ndtr.t.sol b/src/test/Ndtr.t.sol
index 41b15fa..5ea66ce 100644
--- a/src/test/Ndtr.t.sol
+++ b/src/test/Ndtr.t.sol
@@ -206,50 +206,51 @@ contract TestNdtr is Test {
         assertTrue(a >= b, "ndtr-monotonically-decreasing"); // note: not strict?
     }
 
-    /// todo: update this test... it will fail because ndtr is better than the reference, so not equal!
-    function testFuzz_ndtr_reference(int256 x) public {
-        x = bound(x, -int256(2 * 1e27), int256(2 * 1e27));
-        int256 y0 = Ndtr.ndtr(x);
-        x /= 1e9;
-        int256 y1 = Gaussian.cdf(x);
-        console.logInt(y0);
-        console.logInt(y1);
-
-        assertEqPrecision(y0, y1 * 1e9, 1e9, "ndtr-reference");
-    }
-
-    function testFuzz_erfc_ReturnsTwoWhenInputIsTooLow(int256 x) public {
-        vm.assume(x <= Gaussian.ERFC_DOMAIN_LOWER);
-
-        // TODO: Investigate why the error selector is 0x4d2d75b1
-        // instead of 0x35278d12
-        int256 y = Gaussian.erfc(x);
-        assertEq(y, int256(2 ether), "erfc-not-two");
+    function testFuzz_erfc_ReturnsTwoWhenInputIsTooLow(int128 x) public {
+        vm.assume(x < 0);
+        int256 z = -mulfp(x, x);
+        vm.assume(z < -RAY_MAXLOG);
+        int256 y = Ndtr.erfc(x);
+        assertEq(y, int256(RAY_TWO), "erfc-not-two");
     }
 
-    function testFuzz_erfc_ReturnsZeroWhenInputIsTooHigh(int256 x) public {
-        vm.assume(x >= Gaussian.ERFC_DOMAIN_UPPER);
-        int256 y = Gaussian.erfc(x);
+    function testFuzz_erfc_ReturnsZeroWhenInputIsTooHigh(int128 x) public {
+        vm.assume(x > 0);
+        int256 z = -mulfp(x, x);
+        vm.assume(z < -RAY_MAXLOG);
+        int256 y = Ndtr.erfc(x);
         assertEq(y, 0, "erfc-not-zero");
     }
 
     function testFuzz_erfc_NegativeInputIsBounded(int256 x) public {
-        vm.assume(x > -1999999999999999998000000000000000002);
-        vm.assume(x < -0.0000001 ether);
-        int256 y = Gaussian.erfc(x);
-        assertGe(y, 1 ether);
-        assertLe(y, 2 ether);
+        vm.assume(x > -1999999999999999998000000000000000002 * 1e9);
+        vm.assume(x < -0.0000001 ether * 1e9);
+        int256 y = Ndtr.erfc(x);
+        assertGe(y, RAY);
+        assertLe(y, RAY_TWO);
     }
 
     function testFuzz_erfc_PositiveInputIsBounded(int256 x) public {
-        vm.assume(x > 0.0000001 ether);
-        vm.assume(x < 1999999999999999998000000000000000002);
-        int256 y = Gaussian.erfc(x);
+        vm.assume(x > 0.0000001 ether * 1e9);
+        vm.assume(x < 1999999999999999998000000000000000002 * 1e9);
+        int256 y = Ndtr.erfc(x);
         assertGe(y, 0 ether);
-        assertLe(y, 1 ether);
+        assertLe(y, RAY);
     }
 
     function test_erfc_zero_input_returns_one() public {
-        assertEq(Gaussian.erfc(0), 1 ether);
+        assertEq(Ndtr.erfc(0), RAY);
     }
+
+    /// todo: update this test... it will fail because ndtr is better than the reference, so not equal!
+    /* function testFuzz_ndtr_reference(int256 x) public {
+        x = bound(x, -int256(2 * 1e27), int256(2 * 1e27));
+        int256 y0 = Ndtr.ndtr(x);
+        x /= 1e9;
+        int256 y1 = Gaussian.cdf(x);
+        console.logInt(y0);
+        console.logInt(y1);
+
+        assertEqPrecision(y0, y1 * 1e9, 1e9, "ndtr-reference");
+    } */
 }
diff --git a/src/test/RMM01Monotonicity.t.sol b/src/test/RMM01Monotonicity.t.sol
index a2aae65..2dee137 100644
--- a/src/test/RMM01Monotonicity.t.sol
+++ b/src/test/RMM01Monotonicity.t.sol
@@ -5,6 +5,7 @@ import "forge-std/Test.sol";
 import "../Ndtr.sol";
 import "../Ndtri.sol";
 
+/// @dev Non-strict monotonicity at 1E18 precision. Loss of monotonicity beyond 1E18.
 library RMM01 {
     uint256 constant PRECISION = 1e27;
     int256 internal constant YEAR = 31556952;
@@ -31,17 +32,11 @@ library RMM01 {
         } else {
             int256 tauSeconds = divfp(int256(timeRemainingSeconds), YEAR);
             int256 sqrtTau = sqrtfp(tauSeconds); // * 1e18; // sqrt(1e27) gives us a valuew with 1e9 units... so mulfp by 1e18 to get 1e27 units
-            logger.logInt(sqrtTau);
             int256 sqrtTauVol = mulfp(int256(volatilityRay), sqrtTau);
-            logger.logInt(sqrtTauVol);
             int256 phi = int256(PRECISION) - int256(x);
-            logger.logInt(phi);
             int256 inverseCdf = Ndtri.ndtri(phi);
-            logger.logInt(inverseCdf);
             int256 input = inverseCdf - sqrtTauVol;
-            logger.logInt(input);
             int256 cdf = Ndtr.ndtr(input);
-            logger.logInt(cdf);
             int256 product = mulfp(int256(strikePriceRay), cdf);
             y = uint256(product + currentInvariant);
         }
@@ -63,51 +58,126 @@ library RMM01 {
 }
 
 contract TestRMM01Monotonicity is Test {
+    bool DEBUG = false;
+
+    uint256 constant DESIRED_PRECISION = 1e18;
+    uint256 constant DESIRED_PRECISION_SCALAR = 1e9; // Desired precision is 1e18, therefore scalar is 1e27 / 1e18 = 1e9
+
+    /// @dev Asserts `a` is greater than or equal to `b` up to `precision` decimal places.
+    /// @param a The expected larger value.
+    /// @param b The expected smaller value.
+    /// @param precision The number of decimal places to check.
+    /// @param message The message to display if the assertion fails.
+    function assertTruePrecisionNotStrictDecreasing(uint256 a, uint256 b, uint256 precision, string memory message)
+        internal
+    {
+        // Amount to divide by to scale to precision.
+        uint256 scalar = uint256(RAY) / precision;
+
+        // Quantities beyond the radix point at `precision`. These values are truncated before being checked in this assertion.
+        uint256 remainder0 = mulmod(uint256(a), uint256(precision), uint256(RAY));
+        uint256 remainder1 = mulmod(uint256(b), uint256(precision), uint256(RAY));
+
+        if (DEBUG) {
+            console.log("===== RAW AMOUNTS =====");
+            console.log("a", a);
+            console.log("b", b);
+
+            console.log("===== SCALED AMOUNTS =====");
+            console.log("a", a / scalar);
+            console.log("b", b / scalar);
+
+            console.log("===== REMAINDERS =====");
+            console.log("remainder0", remainder0);
+            console.log("remainder1", remainder1);
+        }
+
+        assertTrue(b / scalar <= a / scalar, "b <= a");
+    }
+
+    /// @dev Asserts `b` is greater than or equal to `a` up to `precision` decimal places.
+    /// @param a The expected smaller value.
+    /// @param b The expected larger value.
+    /// @param precision The number of decimal places to check.
+    /// @param message The message to display if the assertion fails.
+    function assertTruePrecisionNotStrictIncreasing(uint256 a, uint256 b, uint256 precision, string memory message)
+        internal
+    {
+        // Amount to divide by to scale to precision.
+        uint256 scalar = uint256(RAY) / precision;
+
+        // Quantities beyond the radix point at `precision`. These values are truncated before being checked in this assertion.
+        uint256 remainder0 = mulmod(uint256(a), uint256(precision), uint256(RAY));
+        uint256 remainder1 = mulmod(uint256(b), uint256(precision), uint256(RAY));
+
+        if (DEBUG) {
+            console.log("===== RAW AMOUNTS =====");
+            console.log("a", a);
+            console.log("b", b);
+
+            console.log("===== SCALED AMOUNTS =====");
+            console.log("a", a / scalar);
+            console.log("b", b / scalar);
+
+            console.log("===== REMAINDERS =====");
+            console.log("remainder0", remainder0);
+            console.log("remainder1", remainder1);
+        }
+
+        assertTrue(b / scalar >= a / scalar, "b >= a");
+    }
+
+    /// ("✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨");
+    /// ("✨✨INVARIANT HAS INCREASING MONOTONICITY✨✨");
+    /// ("✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨");
     /// As x -> 1, y -> 0
-    function testFuzz_rmm_montonically_increasing(uint256 x, uint256 v, uint256 t) public {
-        x = bound(x, 2, 1e27 - 2); // ray
+    function testFuzz_rmm_montonically_increasing(uint256 x1, uint256 strike, uint256 v, uint256 t) public {
+        x1 = bound(x1, DESIRED_PRECISION_SCALAR + 1, 1e27 - 2); // ray
+        strike = bound(strike, DESIRED_PRECISION_SCALAR + 1, 1e27 - 2); // bounds between the lowest value on the desired precision scale and the upper bound of the domain.
         v = bound(v, 1, 1e7); // bps
         t = bound(t, 1, 500 days); // seconds
 
         // scales v up to ray since we bound by a small range in units of bps, since this is a realistic range
         v = v * 1e27 / 1e4; // Divides by BPS units which is 1e4, since 1e4 == 100% == 1.0 == 1e27 ray.
 
-        uint256 strikePriceRay = 1e27; // simple strike price since the domain is large.
-        int256 k = 0;
+        int256 k = 0; // invariant
+        uint256 y1 = RMM01.getY(x1, strike, v, t, k);
+        uint256 y2 = RMM01.getY(x1 + 1, strike, v, t, k);
+        vm.assume(y1 > 0);
 
-        uint256 y = RMM01.getY(x, strikePriceRay, v, t, k);
-        uint256 y2 = RMM01.getY(x + 1, strikePriceRay, v, t, k);
-        vm.assume(y > 0);
-
-        console.log("x", x);
-        console.log("y", y);
-        console.log("x2", x + 1);
+        console.log("x1", x1);
+        console.log("y1", y1);
+        console.log("x2", x1 + 1);
         console.log("y2", y2);
 
-        assertTrue(y <= y2, "y <= y2");
+        // Expect y1 is smaller & y2 is larger, or equal.
+        assertTruePrecisionNotStrictIncreasing(y1, y2, DESIRED_PRECISION, "y2 <= y1");
     }
 
+    /// ("✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨");
+    /// ("✨✨INVARIANT HAS DECREASING MONOTONICITY✨✨");
+    /// ("✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨");
     /// As x -> 0, y -> K
-    function testFuzz_rmm_montonically_decreasing(uint256 x, uint256 v, uint256 t) public {
-        x = bound(x, 2, 1e27 - 2); // ray
+    function testFuzz_rmm_montonically_decreasing(uint256 x1, uint256 strike, uint256 v, uint256 t) public {
+        x1 = bound(x1, DESIRED_PRECISION_SCALAR + 1, 1e27 - 2); // bounds between the lowest value on the desired precision scale and the upper bound of the domain.
+        strike = bound(strike, DESIRED_PRECISION_SCALAR + 1, 1e27 - 2); // bounds between the lowest value on the desired precision scale and the upper bound of the domain.
         v = bound(v, 1, 1e7); // bps
         t = bound(t, 1, 500 days); // seconds
 
         // scales v up to ray since we bound by a small range in units of bps, since this is a realistic range
         v = v * 1e27 / 1e4; // Divides by BPS units which is 1e4, since 1e4 == 100% == 1.0 == 1e27 ray.
 
-        uint256 strikePriceRay = 1e27; // simple strike price since the domain is large.
-        int256 k = 0;
-
-        uint256 y = RMM01.getY(x, strikePriceRay, v, t, k);
-        uint256 y2 = RMM01.getY(x - 1, strikePriceRay, v, t, k);
-        vm.assume(y > 0);
+        int256 k = 0; // invariant
+        uint256 y1 = RMM01.getY(x1, strike, v, t, k);
+        uint256 y2 = RMM01.getY(x1 - 1, strike, v, t, k);
+        vm.assume(y1 > 0);
 
-        console.log("x", x);
-        console.log("y", y);
-        console.log("x2", x - 1);
+        console.log("x1", x1);
+        console.log("y1", y1);
+        console.log("x2", x1 - 1);
         console.log("y2", y2);
 
-        assertTrue(y2 <= y, "y2 <= y");
+        // Expect y1 is larger & y2 is smaller, or equal.
+        assertTruePrecisionNotStrictDecreasing(y1, y2, DESIRED_PRECISION, "y2 <= y1");
     }
 }

From d434e10bf28a5743b8246d80545bbfb9f2a00baf Mon Sep 17 00:00:00 2001
From: alex <alexangel.j@gmail.com>
Date: Sun, 21 May 2023 22:47:09 -0700
Subject: [PATCH 5/9] chore(readme): updates readme cli config guide

---
 README.md | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/README.md b/README.md
index b8301af..bbc4841 100644
--- a/README.md
+++ b/README.md
@@ -2,7 +2,7 @@
 
 > This library is in beta. It's not ready for production.
 
-![](https://visitor-badge.laobi.icu/badge?page_id=solstat) [![](https://github.com/primitivefinance/solstat/actions/workflows/ci.yaml/badge.svg)](https://github.com/primitivefinance/solstat/actions/workflows/ci.yaml) [![](https://dcbadge.vercel.app/api/server/primitive?style=flat)](https://discord.gg/primitive) [![Twitter Badge](https://badgen.net/badge/icon/twitter?icon=twitter&label)](https://twitter.com/primitivefi) 
+![](https://visitor-badge.laobi.icu/badge?page_id=solstat) [![](https://github.com/primitivefinance/solstat/actions/workflows/ci.yaml/badge.svg)](https://github.com/primitivefinance/solstat/actions/workflows/ci.yaml) [![](https://dcbadge.vercel.app/api/server/primitive?style=flat)](https://discord.gg/primitive) [![Twitter Badge](https://badgen.net/badge/icon/twitter?icon=twitter&label)](https://twitter.com/primitivefi)
 
 # SolStat
 
@@ -25,6 +25,7 @@ To compute values using the gaussian.js library, you can use this commmand in th
 yarn cli --cdf {value}
 yarn cli --pdf {value}
 yarn cli --ppf {value}
+yarn cli --ndtr {value}
 ```
 
 ## Irrational Functions

From bf624b058e22ad95f005c7bc4da61fb6772b78d3 Mon Sep 17 00:00:00 2001
From: alex <alexangel.j@gmail.com>
Date: Mon, 22 May 2023 16:29:34 -0700
Subject: [PATCH 6/9] test(ndtr): fix ndtr tests

---
 src/Ndtr.sol                     |   2 +-
 src/test/Ndtr.t.sol              | 194 +++++++++++++++++++------------
 src/test/RMM01Monotonicity.t.sol |  38 +++---
 3 files changed, 138 insertions(+), 96 deletions(-)

diff --git a/src/Ndtr.sol b/src/Ndtr.sol
index dc125b4..8180a7f 100644
--- a/src/Ndtr.sol
+++ b/src/Ndtr.sol
@@ -211,7 +211,7 @@ library Ndtr {
         int256[] memory input1 = copy5(abi.decode(U, (int256[5])));
         // p = T[0] + T[1]*z + ... + T[4]*z^4
         // q = 1 + U[0]*z + ... + U[4]*z^5
-        // y = z * p / q
+        // y = x * p / q
         y = x * polevl(z, input0, 4) / p1evl(z, input1, 5);
 
         return y;
diff --git a/src/test/Ndtr.t.sol b/src/test/Ndtr.t.sol
index 5ea66ce..c213f6a 100644
--- a/src/test/Ndtr.t.sol
+++ b/src/test/Ndtr.t.sol
@@ -3,15 +3,19 @@ pragma solidity ^0.8.4;
 
 import "forge-std/Test.sol";
 import "../Ndtr.sol";
-
 import "../Gaussian.sol";
 
 /// @dev Hardcoded expected values computed manually using `normalcdlower` operation in https://keisan.casio.com/calculator
 contract TestNdtr is Test {
+    uint256 constant MAX_ROUNDING_DELTA = 1;
     int256 constant NDTR_ACCURACY = 1e18;
 
     /// @notice Compares two in256 values up to a precision with a base of RAY.
-    function assertEqPrecision(int256 a, int256 b, int256 precision, string memory message) internal {
+    /// @dev IMPORTANT! Asserts `a` and `b` are within 1 wei up to `precision`.
+    function assertApproxEqPrecision(int256 a, int256 b, int256 precision, string memory message) internal {
+        // Amount to divide by to scale to precision.
+        int256 scalar = RAY / precision;
+
         // Gets the digits passed the precision end point.
         uint256 remainder0 = mulmod(uint256(a), uint256(precision), uint256(RAY));
         uint256 remainder1 = mulmod(uint256(b), uint256(precision), uint256(RAY));
@@ -20,110 +24,117 @@ contract TestNdtr is Test {
         remainder0++;
         remainder1++;
 
+        if (false) {
+            console.log("===== RAW AMOUNTS =====");
+            console.logInt(a);
+            console.logInt(b);
+
+            console.log("===== SCALED AMOUNTS =====");
+            console.logInt(a / scalar);
+            console.logInt(b / scalar);
+
+            console.log("===== REMAINDERS =====");
+            console.log("remainder0", remainder0);
+            console.log("remainder1", remainder1);
+        }
+
         // Converts units to precision.
         a = a * precision / RAY;
         b = b * precision / RAY;
 
-        // Rounds up if remainder is >= 0.5.
-        if (int256(remainder0) >= RAY_HALF) a++;
-        if (int256(remainder1) >= RAY_HALF) b++;
-
-        assertEq(a, b, message);
+        assertApproxEqAbs(a, b, MAX_ROUNDING_DELTA, message);
     }
 
     function test_ndtr() public {
-        assertEqPrecision(Ndtr.ndtr(RAY), int256(0.841344746068542948585232545e27), NDTR_ACCURACY, "ndtr-1");
-        assertEqPrecision(Ndtr.ndtr(RAY * 15 / 10), int256(0.933192798731141934488289354e27), NDTR_ACCURACY, "ndtr-1.5");
-
-        assertEqPrecision(Ndtr.ndtr(RAY * 2), int256(0.977249868051820792947945109e27), NDTR_ACCURACY, "ndtr-2");
-        assertEqPrecision(Ndtr.ndtr(RAY * 10), int256(0.999999999999999999999999999e27), NDTR_ACCURACY, "ndtr-10");
-        assertEqPrecision(
+        assertApproxEqPrecision(Ndtr.ndtr(RAY), int256(0.841344746068542948585232545e27), NDTR_ACCURACY, "ndtr-1");
+        assertApproxEqPrecision(
+            Ndtr.ndtr(RAY * 15 / 10), int256(0.933192798731141934488289354e27), NDTR_ACCURACY, "ndtr-1.5"
+        );
+        assertApproxEqPrecision(Ndtr.ndtr(RAY * 2), int256(0.977249868051820792947945109e27), NDTR_ACCURACY, "ndtr-2");
+        assertApproxEqPrecision(Ndtr.ndtr(RAY * 10), int256(0.999999999999999999999999999e27), NDTR_ACCURACY, "ndtr-10");
+        assertApproxEqPrecision(
             Ndtr.ndtr(RAY / 10000000), int256(0.500000039894228040143796598e27), NDTR_ACCURACY, "ndtr-0.0000001"
         );
-        assertEqPrecision(Ndtr.ndtr(RAY / 100), int256(0.503989356314631603924543868e27), NDTR_ACCURACY, "ndtr-0.01");
+        assertApproxEqPrecision(
+            Ndtr.ndtr(RAY / 100), int256(0.503989356314631603924543868e27), NDTR_ACCURACY, "ndtr-0.01"
+        );
         // negative
-        assertEqPrecision(
+        assertApproxEqPrecision(
             Ndtr.ndtr(-RAY / 10000000), int256(0.499999960105771959856203402e27), NDTR_ACCURACY, "ndtr--0.0000001"
         );
-        assertEqPrecision(Ndtr.ndtr(-RAY / 100), int256(0.496010643685368396075456132e27), NDTR_ACCURACY, "ndtr--0.01");
+        assertApproxEqPrecision(
+            Ndtr.ndtr(-RAY / 100), int256(0.496010643685368396075456132e27), NDTR_ACCURACY, "ndtr--0.01"
+        );
         // 0
-        assertEqPrecision(Ndtr.ndtr(0), int256(0.5e27), NDTR_ACCURACY, "ndtr-0");
+        assertApproxEqPrecision(Ndtr.ndtr(0), int256(0.5e27), NDTR_ACCURACY, "ndtr-0");
         // random values at scale of 1e27
-        assertEqPrecision(
+        assertApproxEqPrecision(
             Ndtr.ndtr(0.123456789e27), int256(0.5491273050781420888711e27), NDTR_ACCURACY, "ndtr-0.123456789"
         );
-        assertEqPrecision(
+        assertApproxEqPrecision(
             Ndtr.ndtr(0.987654321e27), int256(0.83833901356624443490786e27), NDTR_ACCURACY, "ndtr-0.987654321"
         );
-
-        assertEqPrecision(
+        assertApproxEqPrecision(
             Ndtr.ndtr(-0.123456789e27), int256(0.4508726949218579111288e27), NDTR_ACCURACY, "ndtr--0.123456789"
         );
-        assertEqPrecision(
+        assertApproxEqPrecision(
             Ndtr.ndtr(-0.987654321e27), int256(0.16166098643375556509213e27), NDTR_ACCURACY, "ndtr--0.987654321"
         );
     }
 
     /// @dev Tests the edge case where the input is an integer near zero on the ray scale
     function test_ndtr_near_zero() public {
-        assertEqPrecision(Ndtr.ndtr(RAY / NDTR_ACCURACY), int256(0.5e27), NDTR_ACCURACY, "ndtr-1e-18");
-        assertEqPrecision(
+        assertApproxEqPrecision(Ndtr.ndtr(RAY / NDTR_ACCURACY), int256(0.5e27), NDTR_ACCURACY, "ndtr-1e-18");
+        assertApproxEqPrecision(
             Ndtr.ndtr(-RAY / NDTR_ACCURACY), int256(0.499999999999999999999999999e27), NDTR_ACCURACY, "ndtr--1e-18"
         );
-
-        assertEqPrecision(Ndtr.ndtr(RAY / NDTR_ACCURACY + 1), int256(0.5e27), NDTR_ACCURACY, "ndtr-1e-18+1");
-        assertEqPrecision(
+        assertApproxEqPrecision(Ndtr.ndtr(RAY / NDTR_ACCURACY + 1), int256(0.5e27), NDTR_ACCURACY, "ndtr-1e-18+1");
+        assertApproxEqPrecision(
             Ndtr.ndtr(-RAY / NDTR_ACCURACY - 1),
             int256(0.499999999999999999999999999e27),
             NDTR_ACCURACY,
             "ndtr--1e-18-1"
         );
-
-        assertEqPrecision(
+        assertApproxEqPrecision(
             Ndtr.ndtr(RAY / NDTR_ACCURACY * 1e1), int256(0.500000000000000003989422804e27), NDTR_ACCURACY, "ndtr-1e-17"
         );
-        assertEqPrecision(
+        assertApproxEqPrecision(
             Ndtr.ndtr(-RAY / NDTR_ACCURACY * 1e1),
             int256(0.499999999999999996010577195e27),
             NDTR_ACCURACY,
             "ndtr--1e-17"
         );
-
-        assertEqPrecision(
+        assertApproxEqPrecision(
             Ndtr.ndtr(RAY / NDTR_ACCURACY * 1e2), int256(0.50000000000000003989422804e27), NDTR_ACCURACY, "ndtr-1e-16"
         );
-        assertEqPrecision(
+        assertApproxEqPrecision(
             Ndtr.ndtr(-RAY / NDTR_ACCURACY * 1e2),
             int256(0.499999999999999960105771959e27),
             NDTR_ACCURACY,
             "ndtr--1e-16"
         );
-
-        assertEqPrecision(
+        assertApproxEqPrecision(
             Ndtr.ndtr(RAY / NDTR_ACCURACY * 1e3), int256(0.500000000000000398942280401e27), NDTR_ACCURACY, "ndtr-1e-15"
         );
-        assertEqPrecision(
+        assertApproxEqPrecision(
             Ndtr.ndtr(-RAY / NDTR_ACCURACY * 1e3),
             int256(0.499999999999999601057719598e27),
             NDTR_ACCURACY,
             "ndtr--1e-15"
         );
-
-        assertEqPrecision(
+        assertApproxEqPrecision(
             Ndtr.ndtr(RAY / NDTR_ACCURACY * 1e4), int256(0.500000000000003989422804014e27), NDTR_ACCURACY, "ndtr-1e-14"
         );
-        assertEqPrecision(
+        assertApproxEqPrecision(
             Ndtr.ndtr(-RAY / NDTR_ACCURACY * 1e4),
             int256(0.499999999999996010577195985e27),
             NDTR_ACCURACY,
             "ndtr--1e-14"
         );
-
-        assertEqPrecision(
+        assertApproxEqPrecision(
             Ndtr.ndtr(RAY / NDTR_ACCURACY * 1e5), int256(0.500000000000039894228040143e27), NDTR_ACCURACY, "ndtr-1e-13"
         );
-
-        assertEqPrecision(
+        assertApproxEqPrecision(
             Ndtr.ndtr(-RAY / NDTR_ACCURACY * 1e5),
             int256(0.499999999999960105771959856e27),
             NDTR_ACCURACY,
@@ -132,49 +143,50 @@ contract TestNdtr is Test {
     }
 
     function test_ndtr_far_zero() public {
-        assertEqPrecision(Ndtr.ndtr(RAY * 2), int256(0.977249868051820792947945109e27), NDTR_ACCURACY, "ndtr-2");
-        assertEqPrecision(Ndtr.ndtr(RAY * 3), int256(0.998650101968369905473348185e27), NDTR_ACCURACY, "ndtr-3");
-        assertEqPrecision(Ndtr.ndtr(RAY * 4), int256(0.999968328758166880078746229e27), NDTR_ACCURACY, "ndtr-4");
-        assertEqPrecision(Ndtr.ndtr(RAY * 5), int256(0.999999713348428120806088326e27), NDTR_ACCURACY, "ndtr-5");
-        assertEqPrecision(Ndtr.ndtr(RAY * 6), int256(0.999999999013412354962301859e27), NDTR_ACCURACY, "ndtr-6");
-        assertEqPrecision(Ndtr.ndtr(RAY * 7), int256(0.999999999998720187456114164e27), NDTR_ACCURACY, "ndtr-7");
-        assertEqPrecision(Ndtr.ndtr(RAY * 8), int256(0.999999999999999377903942572e27), NDTR_ACCURACY, "ndtr-8");
-        assertEqPrecision(Ndtr.ndtr(RAY * 9), int256(0.999999999999999999887141159e27), NDTR_ACCURACY, "ndtr-9");
-        assertEqPrecision(Ndtr.ndtr(RAY * 10), int256(0.999999999999999999887141159e27), NDTR_ACCURACY, "ndtr-10");
+        assertApproxEqPrecision(Ndtr.ndtr(RAY * 2), int256(0.977249868051820792947945109e27), NDTR_ACCURACY, "ndtr-2");
+        assertApproxEqPrecision(Ndtr.ndtr(RAY * 3), int256(0.998650101968369905473348185e27), NDTR_ACCURACY, "ndtr-3");
+        assertApproxEqPrecision(Ndtr.ndtr(RAY * 4), int256(0.999968328758166880078746229e27), NDTR_ACCURACY, "ndtr-4");
+        assertApproxEqPrecision(Ndtr.ndtr(RAY * 5), int256(0.999999713348428120806088326e27), NDTR_ACCURACY, "ndtr-5");
+        assertApproxEqPrecision(Ndtr.ndtr(RAY * 6), int256(0.999999999013412354962301859e27), NDTR_ACCURACY, "ndtr-6");
+        assertApproxEqPrecision(Ndtr.ndtr(RAY * 7), int256(0.999999999998720187456114164e27), NDTR_ACCURACY, "ndtr-7");
+        assertApproxEqPrecision(Ndtr.ndtr(RAY * 8), int256(0.999999999999999377903942572e27), NDTR_ACCURACY, "ndtr-8");
+        assertApproxEqPrecision(Ndtr.ndtr(RAY * 9), int256(0.999999999999999999887141159e27), NDTR_ACCURACY, "ndtr-9");
+        assertApproxEqPrecision(Ndtr.ndtr(RAY * 10), int256(0.999999999999999999887141159e27), NDTR_ACCURACY, "ndtr-10");
 
         // Negative
-        assertEqPrecision(Ndtr.ndtr(-RAY * 2), int256(0.022750131948179207200282637e27), NDTR_ACCURACY, "ndtr--2");
-        assertEqPrecision(Ndtr.ndtr(-RAY * 3), int256(0.001349898031630094526651814e27), NDTR_ACCURACY, "ndtr--3");
-        assertEqPrecision(Ndtr.ndtr(-RAY * 4), int256(3.167124183311992125377e22), NDTR_ACCURACY, "ndtr--4");
-        assertEqPrecision(Ndtr.ndtr(-RAY * 5), int256(2.86651571879193911673e20), NDTR_ACCURACY, "ndtr--5");
-        assertEqPrecision(Ndtr.ndtr(-RAY * 6), int256(9.8658764503769814e17), NDTR_ACCURACY, "ndtr--6");
-        assertEqPrecision(Ndtr.ndtr(-RAY * 7), int256(1.279812543885835e15), NDTR_ACCURACY, "ndtr--7");
-        assertEqPrecision(Ndtr.ndtr(-RAY * 8), int256(6.22096057427e11), NDTR_ACCURACY, "ndtr--8");
-        assertEqPrecision(Ndtr.ndtr(-RAY * 9), int256(1.1285884e8), NDTR_ACCURACY, "ndtr--9");
-        assertEqPrecision(Ndtr.ndtr(-RAY * 10), int256(7.619e3), NDTR_ACCURACY, "ndtr--10");
+        assertApproxEqPrecision(Ndtr.ndtr(-RAY * 2), int256(0.022750131948179207200282637e27), NDTR_ACCURACY, "ndtr--2");
+        assertApproxEqPrecision(Ndtr.ndtr(-RAY * 3), int256(0.001349898031630094526651814e27), NDTR_ACCURACY, "ndtr--3");
+        assertApproxEqPrecision(Ndtr.ndtr(-RAY * 4), int256(3.167124183311992125377e22), NDTR_ACCURACY, "ndtr--4");
+        assertApproxEqPrecision(Ndtr.ndtr(-RAY * 5), int256(2.86651571879193911673e20), NDTR_ACCURACY, "ndtr--5");
+        assertApproxEqPrecision(Ndtr.ndtr(-RAY * 6), int256(9.8658764503769814e17), NDTR_ACCURACY, "ndtr--6");
+        assertApproxEqPrecision(Ndtr.ndtr(-RAY * 7), int256(1.279812543885835e15), NDTR_ACCURACY, "ndtr--7");
+        assertApproxEqPrecision(Ndtr.ndtr(-RAY * 8), int256(6.22096057427e11), NDTR_ACCURACY, "ndtr--8");
+        assertApproxEqPrecision(Ndtr.ndtr(-RAY * 9), int256(1.1285884e8), NDTR_ACCURACY, "ndtr--9");
+        assertApproxEqPrecision(Ndtr.ndtr(-RAY * 10), int256(7.619e3), NDTR_ACCURACY, "ndtr--10");
     }
 
     /// @dev Expected `b` values computed manually using `normalcdlower` operation in https://keisan.casio.com/calculator
     function test_ndtr_near_half() public {
-        assertEqPrecision(Ndtr.ndtr(RAY / 2), int256(0.69146246127401310363770461e27), NDTR_ACCURACY, "ndtr-0.5");
-        assertEqPrecision(
+        assertApproxEqPrecision(Ndtr.ndtr(RAY / 2), int256(0.69146246127401310363770461e27), NDTR_ACCURACY, "ndtr-0.5");
+        assertApproxEqPrecision(
             Ndtr.ndtr(RAY * 7 / 12), int256(0.720165536400294250547948904e27), NDTR_ACCURACY, "ndtr-0.583333"
         );
-        assertEqPrecision(
+        assertApproxEqPrecision(
             Ndtr.ndtr(RAY * 2 / 3), int256(0.747507462453077086935938175e27), NDTR_ACCURACY, "ndtr-0.666666"
         );
-        assertEqPrecision(
+        assertApproxEqPrecision(
             Ndtr.ndtr(RAY * 5 / 6), int256(0.7976716190363569746314664e27), NDTR_ACCURACY, "ndtr-0.833333"
         );
-
-        assertEqPrecision(Ndtr.ndtr(-RAY / 2), int256(0.308537538725986896362295389e27), NDTR_ACCURACY, "ndtr--0.5");
-        assertEqPrecision(
+        assertApproxEqPrecision(
+            Ndtr.ndtr(-RAY / 2), int256(0.308537538725986896362295389e27), NDTR_ACCURACY, "ndtr--0.5"
+        );
+        assertApproxEqPrecision(
             Ndtr.ndtr(-RAY * 7 / 12), int256(0.279834463599705749452051095e27), NDTR_ACCURACY, "ndtr--0.583333"
         );
-        assertEqPrecision(
+        assertApproxEqPrecision(
             Ndtr.ndtr(-RAY * 2 / 3), int256(0.252492537546922913064061824e27), NDTR_ACCURACY, "ndtr--0.666666"
         );
-        assertEqPrecision(
+        assertApproxEqPrecision(
             Ndtr.ndtr(-RAY * 5 / 6), int256(0.202328380963643025368533599e27), NDTR_ACCURACY, "ndtr--0.833333"
         );
     }
@@ -183,6 +195,7 @@ contract TestNdtr is Test {
     int256 constant LOWER_NDTR_BOUND = -10 * RAY;
 
     /// note: fails strictly at values above abs(1), i.e. a == b
+    /// @dev As x -> inf, y -> 1
     function testFuzz_ndtr_monotonically_increasing(int256 x) public {
         x = bound(x, LOWER_NDTR_BOUND, UPPER_NDTR_BOUND);
 
@@ -195,6 +208,7 @@ contract TestNdtr is Test {
     }
 
     /// note: fails strictly at values above abs(1), i.e. a == b
+    /// @dev As x -> -inf, y -> 0
     function testFuzz_ndtr_monotonically_decreasing(int256 x) public {
         x = bound(x, LOWER_NDTR_BOUND, UPPER_NDTR_BOUND);
 
@@ -206,6 +220,34 @@ contract TestNdtr is Test {
         assertTrue(a >= b, "ndtr-monotonically-decreasing"); // note: not strict?
     }
 
+    /// note: fails strictly at values above abs(1), i.e. a == b
+    /// @dev As x -> -2, y -> 2
+    function testFuzz_erfc_monotonically_increasing(int256 x) public {
+        x = bound(x, LOWER_NDTR_BOUND, UPPER_NDTR_BOUND);
+
+        int256 a = Ndtr.erfc(x);
+        int256 b = Ndtr.erfc(x - 1);
+        console.logInt(a);
+        console.logInt(b);
+
+        assertTrue(a <= b, "erfc-monotonically-increasing"); // note: not strict?
+    }
+
+    /// note: fails strictly at values above abs(1), i.e. a == b
+    /// @dev As x -> 2, y -> 0
+    function testFuzz_erfc_monotonically_decreasing(int256 x) public {
+        vm.assume(x != 0);
+        x = bound(x, LOWER_NDTR_BOUND, UPPER_NDTR_BOUND);
+        x = absolute(x); // erfc is decreasing for positive values
+
+        int256 a = Ndtr.erfc(x);
+        int256 b = Ndtr.erfc(x + 1);
+        console.logInt(a);
+        console.logInt(b);
+
+        assertTrue(a >= b, "erfc-monotonically-decreasing"); // note: not strict?
+    }
+
     function testFuzz_erfc_ReturnsTwoWhenInputIsTooLow(int128 x) public {
         vm.assume(x < 0);
         int256 z = -mulfp(x, x);
@@ -222,17 +264,19 @@ contract TestNdtr is Test {
         assertEq(y, 0, "erfc-not-zero");
     }
 
+    /// @dev erfc domain is -2 < 0 < 2, check between (-2, 0)
     function testFuzz_erfc_NegativeInputIsBounded(int256 x) public {
-        vm.assume(x > -1999999999999999998000000000000000002 * 1e9);
-        vm.assume(x < -0.0000001 ether * 1e9);
+        vm.assume(x > -RAY_TWO); //
+        vm.assume(x < -1e9); // smallest unit of precision closest to zero.
         int256 y = Ndtr.erfc(x);
         assertGe(y, RAY);
         assertLe(y, RAY_TWO);
     }
 
+    /// @dev erfc domain is -2 < 0 < 2, check between (0, 2)
     function testFuzz_erfc_PositiveInputIsBounded(int256 x) public {
-        vm.assume(x > 0.0000001 ether * 1e9);
-        vm.assume(x < 1999999999999999998000000000000000002 * 1e9);
+        vm.assume(x > 1e9);
+        vm.assume(x < RAY_TWO);
         int256 y = Ndtr.erfc(x);
         assertGe(y, 0 ether);
         assertLe(y, RAY);
@@ -251,6 +295,6 @@ contract TestNdtr is Test {
         console.logInt(y0);
         console.logInt(y1);
 
-        assertEqPrecision(y0, y1 * 1e9, 1e9, "ndtr-reference");
+        assertApproxEqPrecision(y0, y1 * 1e9, 1e9, "ndtr-reference");
     } */
 }
diff --git a/src/test/RMM01Monotonicity.t.sol b/src/test/RMM01Monotonicity.t.sol
index 2dee137..e4ce7f5 100644
--- a/src/test/RMM01Monotonicity.t.sol
+++ b/src/test/RMM01Monotonicity.t.sol
@@ -68,9 +68,7 @@ contract TestRMM01Monotonicity is Test {
     /// @param b The expected smaller value.
     /// @param precision The number of decimal places to check.
     /// @param message The message to display if the assertion fails.
-    function assertTruePrecisionNotStrictDecreasing(uint256 a, uint256 b, uint256 precision, string memory message)
-        internal
-    {
+    function assertGtePrecisionNotStrict(uint256 a, uint256 b, uint256 precision, string memory message) internal {
         // Amount to divide by to scale to precision.
         uint256 scalar = uint256(RAY) / precision;
 
@@ -92,7 +90,7 @@ contract TestRMM01Monotonicity is Test {
             console.log("remainder1", remainder1);
         }
 
-        assertTrue(b / scalar <= a / scalar, "b <= a");
+        assertTrue(a / scalar >= b / scalar, "a >= b");
     }
 
     /// @dev Asserts `b` is greater than or equal to `a` up to `precision` decimal places.
@@ -100,9 +98,7 @@ contract TestRMM01Monotonicity is Test {
     /// @param b The expected larger value.
     /// @param precision The number of decimal places to check.
     /// @param message The message to display if the assertion fails.
-    function assertTruePrecisionNotStrictIncreasing(uint256 a, uint256 b, uint256 precision, string memory message)
-        internal
-    {
+    function assertLtePrecisionNotStrict(uint256 a, uint256 b, uint256 precision, string memory message) internal {
         // Amount to divide by to scale to precision.
         uint256 scalar = uint256(RAY) / precision;
 
@@ -124,14 +120,14 @@ contract TestRMM01Monotonicity is Test {
             console.log("remainder1", remainder1);
         }
 
-        assertTrue(b / scalar >= a / scalar, "b >= a");
+        assertTrue(a / scalar <= b / scalar, "a <= b");
     }
 
     /// ("✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨");
-    /// ("✨✨INVARIANT HAS INCREASING MONOTONICITY✨✨");
+    /// ("✨✨INVARIANT HAS DECREASING MONOTONICITY✨✨");
     /// ("✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨");
     /// As x -> 1, y -> 0
-    function testFuzz_rmm_montonically_increasing(uint256 x1, uint256 strike, uint256 v, uint256 t) public {
+    function testFuzz_rmm_montonically_decreasing(uint256 x1, uint256 strike, uint256 v, uint256 t) public {
         x1 = bound(x1, DESIRED_PRECISION_SCALAR + 1, 1e27 - 2); // ray
         strike = bound(strike, DESIRED_PRECISION_SCALAR + 1, 1e27 - 2); // bounds between the lowest value on the desired precision scale and the upper bound of the domain.
         v = bound(v, 1, 1e7); // bps
@@ -141,8 +137,8 @@ contract TestRMM01Monotonicity is Test {
         v = v * 1e27 / 1e4; // Divides by BPS units which is 1e4, since 1e4 == 100% == 1.0 == 1e27 ray.
 
         int256 k = 0; // invariant
-        uint256 y1 = RMM01.getY(x1, strike, v, t, k);
-        uint256 y2 = RMM01.getY(x1 + 1, strike, v, t, k);
+        uint256 y1 = RMM01.getY(x1, strike, v, t, k); // x smaller, y larger
+        uint256 y2 = RMM01.getY(x1 + 1, strike, v, t, k); // x larger, y smaller
         vm.assume(y1 > 0);
 
         console.log("x1", x1);
@@ -150,15 +146,16 @@ contract TestRMM01Monotonicity is Test {
         console.log("x2", x1 + 1);
         console.log("y2", y2);
 
-        // Expect y1 is smaller & y2 is larger, or equal.
-        assertTruePrecisionNotStrictIncreasing(y1, y2, DESIRED_PRECISION, "y2 <= y1");
+        // As x gets larger, expect y to get smaller.
+        // Check that y2 is smaller than y1.
+        assertGtePrecisionNotStrict(y1, y2, DESIRED_PRECISION, "y1 >= y2");
     }
 
     /// ("✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨");
-    /// ("✨✨INVARIANT HAS DECREASING MONOTONICITY✨✨");
+    /// ("✨✨INVARIANT HAS INCREASING MONOTONICITY✨✨");
     /// ("✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨");
     /// As x -> 0, y -> K
-    function testFuzz_rmm_montonically_decreasing(uint256 x1, uint256 strike, uint256 v, uint256 t) public {
+    function testFuzz_rmm_montonically_increasing(uint256 x1, uint256 strike, uint256 v, uint256 t) public {
         x1 = bound(x1, DESIRED_PRECISION_SCALAR + 1, 1e27 - 2); // bounds between the lowest value on the desired precision scale and the upper bound of the domain.
         strike = bound(strike, DESIRED_PRECISION_SCALAR + 1, 1e27 - 2); // bounds between the lowest value on the desired precision scale and the upper bound of the domain.
         v = bound(v, 1, 1e7); // bps
@@ -168,8 +165,8 @@ contract TestRMM01Monotonicity is Test {
         v = v * 1e27 / 1e4; // Divides by BPS units which is 1e4, since 1e4 == 100% == 1.0 == 1e27 ray.
 
         int256 k = 0; // invariant
-        uint256 y1 = RMM01.getY(x1, strike, v, t, k);
-        uint256 y2 = RMM01.getY(x1 - 1, strike, v, t, k);
+        uint256 y1 = RMM01.getY(x1, strike, v, t, k); // x is larger, y is smaller
+        uint256 y2 = RMM01.getY(x1 - 1, strike, v, t, k); // x is smaller, y is larger
         vm.assume(y1 > 0);
 
         console.log("x1", x1);
@@ -177,7 +174,8 @@ contract TestRMM01Monotonicity is Test {
         console.log("x2", x1 - 1);
         console.log("y2", y2);
 
-        // Expect y1 is larger & y2 is smaller, or equal.
-        assertTruePrecisionNotStrictDecreasing(y1, y2, DESIRED_PRECISION, "y2 <= y1");
+        // As x gets smaller, expect y to get larger.
+        // Check that y1 is smaller than y2.
+        assertLtePrecisionNotStrict(y1, y2, DESIRED_PRECISION, "y1 <= y2");
     }
 }

From 8836004719dac7db4e964255acb439a628e1efc6 Mon Sep 17 00:00:00 2001
From: alex <alexangel.j@gmail.com>
Date: Mon, 22 May 2023 17:05:04 -0700
Subject: [PATCH 7/9] test(ndtri): updates ndtri tests and constant

---
 src/Ndtri.sol        |  10 ++-
 src/test/Ndtr.t.sol  |  17 +----
 src/test/Ndtri.t.sol | 150 +++++++++++++++++++++----------------------
 3 files changed, 82 insertions(+), 95 deletions(-)

diff --git a/src/Ndtri.sol b/src/Ndtri.sol
index db213ba..191e81f 100644
--- a/src/Ndtri.sol
+++ b/src/Ndtri.sol
@@ -20,6 +20,12 @@ library Ndtri {
 
     error MaxNumError();
 
+    ////////////////
+    // Constants  //
+    ////////////////
+
+    int256 constant RAY_EXPNEG2 = 0.135335283236612691893999494e27;
+
     //////////////////////////////////////////////
     // Approximation for 0 <= |y - 0.5| <= 3/8  //
     //////////////////////////////////////////////
@@ -140,13 +146,13 @@ library Ndtri {
         y = y0;
 
         // Use a different approximation if y > exp(-2).
-        if (y > RAY_ONE - 0.13533528323661269189e27) {
+        if (y > RAY_ONE - RAY_EXPNEG2) {
             // 0.135... = exp(-2)
             y = RAY_ONE - y;
             code = 0;
         }
 
-        if (y > 0.13533528323661269189e27) {
+        if (y > RAY_EXPNEG2) {
             int256[] memory P0_ARRAY = copy5(abi.decode(P0, (int256[5])));
             int256[] memory Q0_ARRAY = copy8(abi.decode(Q0, (int256[8])));
             y = y - RAY_HALF;
diff --git a/src/test/Ndtr.t.sol b/src/test/Ndtr.t.sol
index c213f6a..ca918da 100644
--- a/src/test/Ndtr.t.sol
+++ b/src/test/Ndtr.t.sol
@@ -20,10 +20,7 @@ contract TestNdtr is Test {
         uint256 remainder0 = mulmod(uint256(a), uint256(precision), uint256(RAY));
         uint256 remainder1 = mulmod(uint256(b), uint256(precision), uint256(RAY));
 
-        // Add one to the remainder to round up, in the case it is 99...99.
-        remainder0++;
-        remainder1++;
-
+        // For debugging...
         if (false) {
             console.log("===== RAW AMOUNTS =====");
             console.logInt(a);
@@ -285,16 +282,4 @@ contract TestNdtr is Test {
     function test_erfc_zero_input_returns_one() public {
         assertEq(Ndtr.erfc(0), RAY);
     }
-
-    /// todo: update this test... it will fail because ndtr is better than the reference, so not equal!
-    /* function testFuzz_ndtr_reference(int256 x) public {
-        x = bound(x, -int256(2 * 1e27), int256(2 * 1e27));
-        int256 y0 = Ndtr.ndtr(x);
-        x /= 1e9;
-        int256 y1 = Gaussian.cdf(x);
-        console.logInt(y0);
-        console.logInt(y1);
-
-        assertApproxEqPrecision(y0, y1 * 1e9, 1e9, "ndtr-reference");
-    } */
 }
diff --git a/src/test/Ndtri.t.sol b/src/test/Ndtri.t.sol
index 08e0f55..89ad6c0 100644
--- a/src/test/Ndtri.t.sol
+++ b/src/test/Ndtri.t.sol
@@ -8,44 +8,60 @@ import "../Gaussian.sol";
 /// @dev Hardcoded expected values computed manually using `normalcdlower` operation in https://keisan.casio.com/calculator
 contract TestNdtri is Test {
     int256 constant NDTRI_ACCURACY = 1e18;
+    uint256 constant MAX_ROUNDING_DELTA = 1;
 
     /// @notice Compares two in256 values up to a precision with a base of RAY.
-    function assertEqPrecision(int256 a, int256 b, int256 precision, string memory message) internal {
+    /// @dev IMPORTANT! Asserts `a` and `b` are within 1 wei up to `precision`.
+    function assertApproxEqPrecision(int256 a, int256 b, int256 precision, string memory message) internal {
+        // Amount to divide by to scale to precision.
+        int256 scalar = RAY / precision;
+
         // Gets the digits passed the precision end point.
         uint256 remainder0 = mulmod(uint256(a), uint256(precision), uint256(RAY));
         uint256 remainder1 = mulmod(uint256(b), uint256(precision), uint256(RAY));
 
-        // Add one to the remainder to round up, in the case it is 99...99.
-        remainder0++;
-        remainder1++;
+        // For debugging...
+        if (false) {
+            console.log("===== RAW AMOUNTS =====");
+            console.logInt(a);
+            console.logInt(b);
+
+            console.log("===== SCALED AMOUNTS =====");
+            console.logInt(a / scalar);
+            console.logInt(b / scalar);
+
+            console.log("===== REMAINDERS =====");
+            console.log("remainder0", remainder0);
+            console.log("remainder1", remainder1);
+        }
 
         // Converts units to precision.
         a = a * precision / RAY;
         b = b * precision / RAY;
 
-        // Rounds up if remainder is >= 0.5.
-        //if (int256(remainder0) >= RAY_HALF) a++;
-        //if (int256(remainder1) >= RAY_HALF) b++;
-
-        assertEq(a, b, message);
+        assertApproxEqAbs(a, b, MAX_ROUNDING_DELTA, message);
     }
 
-    /// @dev Tests various inputs on the domain [0, 1], where input is a probability between 0 and 1.
+    /// @dev Tests various inputs on the domain (0, 1), where input is a probability between 0 and 1.
     function test_ndtri() public {
-        assertEqPrecision(Ndtri.ndtri(RAY / 2), int256(0), NDTRI_ACCURACY, "ndtri-0.5");
-        assertEqPrecision(Ndtri.ndtri(RAY / 4), int256(-0.674489750196081743202227014e27), NDTRI_ACCURACY, "ndtri-0.25");
-        assertEqPrecision(Ndtri.ndtri(RAY / 8), int256(-1.15034938037600817829676531e27), NDTRI_ACCURACY, "ndtri-0.125");
-        assertEqPrecision(
+        assertApproxEqPrecision(Ndtri.ndtri(RAY / 2), int256(0), NDTRI_ACCURACY, "ndtri-0.5");
+        assertApproxEqPrecision(
+            Ndtri.ndtri(RAY / 4), int256(-0.674489750196081743202227014e27), NDTRI_ACCURACY, "ndtri-0.25"
+        );
+        assertApproxEqPrecision(
+            Ndtri.ndtri(RAY / 8), int256(-1.15034938037600817829676531e27), NDTRI_ACCURACY, "ndtri-0.125"
+        );
+        assertApproxEqPrecision(
             Ndtri.ndtri(RAY / 16), int256(-1.53412054435254631170839905e27), NDTRI_ACCURACY, "ndtri-0.0625"
         );
-        assertEqPrecision(
+        assertApproxEqPrecision(
             Ndtri.ndtri(RAY * 77 / 93), int256(0.946122671364591906727098902e27), NDTRI_ACCURACY, "ndtri-0.827956"
         );
-        assertEqPrecision(
+        assertApproxEqPrecision(
             Ndtri.ndtri(RAY * 11 / 85), int256(-1.1291761577077979287868872e27), NDTRI_ACCURACY, "ndtri-0.129412"
         );
 
-        assertEqPrecision(
+        assertApproxEqPrecision(
             Ndtri.ndtri(729329433526774616220000000),
             int256(0.610786196335546300516717471e27),
             NDTRI_ACCURACY,
@@ -53,78 +69,86 @@ contract TestNdtri is Test {
         );
     }
 
-    /// Starting at 1e27, decrements by atleast RAY / NDTRI_ACCURACY = 1e9, since thats the lowest unit
+    /// @dev Starting at 1e27, decrements by atleast RAY / NDTRI_ACCURACY = 1e9, since thats the lowest unit
     /// that has precision.
-    function test_ndtri_near_one() public {
+    /// i.e. x = 1e27 - 1e9 = 0.999999999999999999000000000
+    /// @notice Upper bound is: 1 - exp(-2) < x < 1
+    function test_ndtri_upper_bound() public {
         int256 decrement = RAY / NDTRI_ACCURACY;
-        console.logInt(RAY - decrement);
-        assertEqPrecision(
+        assertApproxEqPrecision(
             Ndtri.ndtri(RAY - decrement),
             int256(8.75729034878231506388112862e27),
             NDTRI_ACCURACY,
             "ndtri-0.999999999999999999000000000"
         );
-        assertEqPrecision(
+        assertApproxEqPrecision(
             Ndtri.ndtri(RAY - decrement * 1e1), int256(8.49379322410959807444471881e27), NDTRI_ACCURACY, "ndtri-1-1e-17"
         );
-        assertEqPrecision(
+        assertApproxEqPrecision(
             Ndtri.ndtri(RAY - decrement * 1e2), int256(8.22208221613043561267585878e27), NDTRI_ACCURACY, "ndtri-1-1e-16"
         );
-        assertEqPrecision(
+        assertApproxEqPrecision(
             Ndtri.ndtri(RAY - decrement * 1e3), int256(7.94134532617099678096674357e27), NDTRI_ACCURACY, "ndtri-1-1e-15"
         );
-        assertEqPrecision(
+        assertApproxEqPrecision(
             Ndtri.ndtri(RAY - decrement * 1e4), int256(7.65062809293526881641896581e27), NDTRI_ACCURACY, "ndtri-1-1e-14"
         );
-        assertEqPrecision(
+        assertApproxEqPrecision(
             Ndtri.ndtri(RAY - decrement * 1e5), int256(7.34879610280067751753910726e27), NDTRI_ACCURACY, "ndtri-1-1e-13"
         );
     }
 
     /// @dev Tests the edge case where the input is an integer near zero on the ray scale
-    function test_ndtri_near_zero() public {
-        assertEqPrecision(
+    /// i.e. x = 1e9 = 0.000000000000000001000000000
+    /// @notice Lower bound is: 0 < x < exp(-2)
+    function test_ndtri_lower_bound() public {
+        console.logInt(Gaussian.ppf(1e1));
+        assertApproxEqPrecision(
             Ndtri.ndtri(RAY / NDTRI_ACCURACY + 1), -int256(8.7572903487823e27), NDTRI_ACCURACY, "ndtri-1e-18"
         );
-        assertEqPrecision(
+        assertApproxEqPrecision(
             Ndtri.ndtri(RAY / NDTRI_ACCURACY * 1e1), -int256(8.493793224109598e27), NDTRI_ACCURACY, "ndtri-1e-17"
         );
-        assertEqPrecision(
+        assertApproxEqPrecision(
             Ndtri.ndtri(RAY / NDTRI_ACCURACY * 1e2), -int256(8.2220822161304356e27), NDTRI_ACCURACY, "ndtri-1e-16"
         );
-        assertEqPrecision(
+        assertApproxEqPrecision(
             Ndtri.ndtri(RAY / NDTRI_ACCURACY * 1e3), -int256(7.9413453261709968e27), NDTRI_ACCURACY, "ndtri-1e-15"
         );
-        assertEqPrecision(
+        assertApproxEqPrecision(
             Ndtri.ndtri(RAY / NDTRI_ACCURACY * 1e4), -int256(7.65062809293526882e27), NDTRI_ACCURACY, "ndtri-1e-14"
         );
-        assertEqPrecision(
+        assertApproxEqPrecision(
             Ndtri.ndtri(RAY / NDTRI_ACCURACY * 1e5), -int256(7.348796102800677518e27), NDTRI_ACCURACY, "ndtri-1e-13"
         );
     }
 
-    function test_ndtri_far_zero() public {}
-
     /// @dev Expected `b` values computed manually using `normalcdlower` operation in https://keisan.casio.com/calculator
-    function test_ndtri_near_half() public {
-        assertEqPrecision(
-            Ndtri.ndtri(RAY * 3 / 5), int256(0.253347103135799798798196181e27), NDTRI_ACCURACY, "ndtri-0.666666"
+    /// i.e. x = 0.5 = 0.500000000000000000000000000
+    /// @notice Middle bound is: exp(-2) < x < 1 - exp(-2)
+    function test_ndtri_middle_bound() public {
+        assertApproxEqPrecision(
+            Ndtri.ndtri(RAY * 1 / 4), int256(-0.674489750196081743202227014e27), NDTRI_ACCURACY, "ndtri-0.25"
+        );
+        assertApproxEqPrecision(
+            Ndtri.ndtri(RAY * 7 / 12), int256(0.210428394247924723324540643e27), NDTRI_ACCURACY, "ndtri-0.583333..."
         );
-        assertEqPrecision(
-            Ndtri.ndtri(RAY * 7 / 12), int256(0.210428394247924723324540643e27), NDTRI_ACCURACY, "ndtri-0.583333"
+        assertApproxEqPrecision(
+            Ndtri.ndtri(RAY * 3 / 5), int256(0.253347103135799798798196181e27), NDTRI_ACCURACY, "ndtri-0.6"
         );
-        assertEqPrecision(
-            Ndtri.ndtri(RAY * 2 / 3), int256(0.43072729929545749020594039e27), NDTRI_ACCURACY, "ndtri-0.666666"
+        assertApproxEqPrecision(
+            Ndtri.ndtri(RAY * 2 / 3), int256(0.43072729929545749020594039e27), NDTRI_ACCURACY, "ndtri-0.666666..."
         );
-        assertEqPrecision(
-            Ndtri.ndtri(RAY * 5 / 6), int256(0.96742156610170103955040122e27), NDTRI_ACCURACY, "ndtri-0.833333"
+        assertApproxEqPrecision(
+            Ndtri.ndtri(RAY * 3 / 4), int256(0.674489750196081743202227014e27), NDTRI_ACCURACY, "ndtri-0.75"
         );
     }
 
     int256 constant UPPER_NDTRI_BOUND = RAY - 1 - 1 - 1; // Bound is RAY - 1, but we want to increase by 1 so need at least 2
     int256 constant LOWER_NDTRI_BOUND = 1 + 1 + 1; // Bound is 1, but we want to decrease by 1, so need at least 2
 
-    /// note: strict passes... does it pass with 10k fuzz runs though?
+    /// note: not strict!
+    /// @dev As x -> 1, ndtri(x) -> inf
     function testFuzz_ndtri_monotonically_increasing(int256 x) public {
         x = bound(x, LOWER_NDTRI_BOUND, UPPER_NDTRI_BOUND);
 
@@ -133,10 +157,11 @@ contract TestNdtri is Test {
         console.logInt(a);
         console.logInt(b);
 
-        assertTrue(a < b, "ndtri-monotonically-increasing");
+        assertTrue(a <= b, "ndtri-monotonically-increasing");
     }
 
-    /// note: strict passes... does it pass with 10k fuzz runs though?
+    /// note: not strict!
+    /// @dev As x -> 0, ndtri(x) -> -inf
     function testFuzz_ndtri_monotonically_decreasing(int256 x) public {
         x = bound(x, LOWER_NDTRI_BOUND, UPPER_NDTRI_BOUND);
 
@@ -145,7 +170,7 @@ contract TestNdtri is Test {
         console.logInt(a);
         console.logInt(b);
 
-        assertTrue(a > b, "ndtri-monotonically-decreasing");
+        assertTrue(a >= b, "ndtri-monotonically-decreasing");
     }
 
     function test_ndtri_equals_one_reverts() public {
@@ -167,33 +192,4 @@ contract TestNdtri is Test {
         vm.expectRevert(Ndtri.MaxNumError.selector);
         Ndtri.ndtri(0 - 1);
     }
-
-    function test_ndtri_lower() public {
-        assertEqPrecision(Ndtri.ndtri(0.1e27), int256(-1.28155156554460046696510332e27), NDTRI_ACCURACY, "ndtri-0.1");
-    }
-
-    /* /// todo: update this test... it will fail because ndtri is better than the reference, so not equal!
-    function testFuzz_ndtri_reference(int256 x) public {
-        x = bound(x, 2, 1e27 - 2);
-        int256 y0 = Ndtri.ndtri(x);
-        x /= 1e9;
-        int256 y1 = Gaussian.ppf(x);
-        console.logInt(y0);
-        console.logInt(y1);
-
-        assertEqPrecision(y0, y1 * 1e9, 1e9, "ndtri-reference");
-    }
-
-    /// todo: update this test... it will fail because ndtri is better than the reference, so not equal!
-    function testFuzz_ndtri_refrence_below_exp2(int256 x) public {
-        x = bound(x, 0.13533528323661269189e27 + 1, 1e27 - 0.13533528323661269189e27);
-
-        int256 y0 = Ndtri.ndtri(x);
-        x /= 1e9;
-        int256 y1 = Gaussian.ppf(x);
-        console.logInt(y0);
-        console.logInt(y1);
-
-        assertEqPrecision(y0, y1 * 1e9, 1e9, "ndtri-reference-below-exp2");
-    } */
 }

From 47bc8ac3704d62f69c26324a9214117b19654c48 Mon Sep 17 00:00:00 2001
From: alex <alexangel.j@gmail.com>
Date: Mon, 22 May 2023 18:22:30 -0700
Subject: [PATCH 8/9] test(input): adds tests for input precision

---
 src/test/Ndtr.t.sol              | 16 +++++++++
 src/test/RMM01Monotonicity.t.sol | 59 ++++++++++++++++++++++++++++++--
 src/test/RayMathLib.t.sol        |  4 ++-
 3 files changed, 76 insertions(+), 3 deletions(-)

diff --git a/src/test/Ndtr.t.sol b/src/test/Ndtr.t.sol
index ca918da..4ec48ec 100644
--- a/src/test/Ndtr.t.sol
+++ b/src/test/Ndtr.t.sol
@@ -282,4 +282,20 @@ contract TestNdtr is Test {
     function test_erfc_zero_input_returns_one() public {
         assertEq(Ndtr.erfc(0), RAY);
     }
+
+    function testFuzz_ndtr_input_precision(int256 x) public {
+        x = bound(x, -10e27, 10e27); // Bound between [-10, 10] in units of 1E27.
+        vm.assume(absolute(x) > 1e9); // Since we are dividing by 1E9... need to have at least 1 in the input, so x has to be gte 1E9.
+
+        // Assume we want to use `ndtr()` but have an input with WAD units (1E18).
+        // The input can be scaled to match the units of 1E27.
+        // But do we lose precision in the output?
+        int256 outputWith1E27Input = Ndtr.ndtr(x);
+        x = x / 1e9; // Scale down to 1E18.
+        x = x * 1e9; // Scale back up to 1E27, losing all precision past the 1E18 radix point.
+        int256 outputWith1E18Input = Ndtr.ndtr(x);
+
+        // Asserts that the outputs scaled to 1E18 are equal up to 1E18 precision, regardless of the input precision.
+        assertApproxEqPrecision(outputWith1E27Input, outputWith1E18Input, 1e18, "ndtr-input-precision");
+    }
 }
diff --git a/src/test/RMM01Monotonicity.t.sol b/src/test/RMM01Monotonicity.t.sol
index e4ce7f5..4709f05 100644
--- a/src/test/RMM01Monotonicity.t.sol
+++ b/src/test/RMM01Monotonicity.t.sol
@@ -60,9 +60,42 @@ library RMM01 {
 contract TestRMM01Monotonicity is Test {
     bool DEBUG = false;
 
+    uint256 constant MAX_ROUNDING_DELTA = 1;
     uint256 constant DESIRED_PRECISION = 1e18;
     uint256 constant DESIRED_PRECISION_SCALAR = 1e9; // Desired precision is 1e18, therefore scalar is 1e27 / 1e18 = 1e9
 
+    /// @notice Compares two in256 values up to a precision with a base of RAY.
+    /// @dev IMPORTANT! Asserts `a` and `b` are within 1 wei up to `precision`.
+    function assertApproxEqPrecision(uint256 a, uint256 b, uint256 precision, string memory message) internal {
+        // Amount to divide by to scale to precision.
+        uint256 scalar = uint256(RAY) / precision;
+
+        // Gets the digits passed the precision end point.
+        uint256 remainder0 = mulmod(uint256(a), uint256(precision), uint256(RAY));
+        uint256 remainder1 = mulmod(uint256(b), uint256(precision), uint256(RAY));
+
+        // For debugging...
+        if (false) {
+            console.log("===== RAW AMOUNTS =====");
+            console.log("a", a);
+            console.log("b", b);
+
+            console.log("===== SCALED AMOUNTS =====");
+            console.log("a / scalar", a / scalar);
+            console.log("b / scalar", b / scalar);
+
+            console.log("===== REMAINDERS =====");
+            console.log("remainder0", remainder0);
+            console.log("remainder1", remainder1);
+        }
+
+        // Converts units to precision.
+        a = a * precision / uint256(RAY);
+        b = b * precision / uint256(RAY);
+
+        assertApproxEqAbs(a, b, MAX_ROUNDING_DELTA, message);
+    }
+
     /// @dev Asserts `a` is greater than or equal to `b` up to `precision` decimal places.
     /// @param a The expected larger value.
     /// @param b The expected smaller value.
@@ -127,7 +160,7 @@ contract TestRMM01Monotonicity is Test {
     /// ("✨✨INVARIANT HAS DECREASING MONOTONICITY✨✨");
     /// ("✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨");
     /// As x -> 1, y -> 0
-    function testFuzz_rmm_montonically_decreasing(uint256 x1, uint256 strike, uint256 v, uint256 t) public {
+    function testFuzz_getY_montonically_decreasing(uint256 x1, uint256 strike, uint256 v, uint256 t) public {
         x1 = bound(x1, DESIRED_PRECISION_SCALAR + 1, 1e27 - 2); // ray
         strike = bound(strike, DESIRED_PRECISION_SCALAR + 1, 1e27 - 2); // bounds between the lowest value on the desired precision scale and the upper bound of the domain.
         v = bound(v, 1, 1e7); // bps
@@ -155,7 +188,7 @@ contract TestRMM01Monotonicity is Test {
     /// ("✨✨INVARIANT HAS INCREASING MONOTONICITY✨✨");
     /// ("✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨");
     /// As x -> 0, y -> K
-    function testFuzz_rmm_montonically_increasing(uint256 x1, uint256 strike, uint256 v, uint256 t) public {
+    function testFuzz_getY_montonically_increasing(uint256 x1, uint256 strike, uint256 v, uint256 t) public {
         x1 = bound(x1, DESIRED_PRECISION_SCALAR + 1, 1e27 - 2); // bounds between the lowest value on the desired precision scale and the upper bound of the domain.
         strike = bound(strike, DESIRED_PRECISION_SCALAR + 1, 1e27 - 2); // bounds between the lowest value on the desired precision scale and the upper bound of the domain.
         v = bound(v, 1, 1e7); // bps
@@ -178,4 +211,26 @@ contract TestRMM01Monotonicity is Test {
         // Check that y1 is smaller than y2.
         assertLtePrecisionNotStrict(y1, y2, DESIRED_PRECISION, "y1 <= y2");
     }
+
+    function testFuzz_getY_input_precision(uint256 x1, uint256 strike, uint256 v, uint256 t) public {
+        x1 = bound(x1, DESIRED_PRECISION_SCALAR + 1, 1e27 - 2); // bounds between the lowest value on the desired precision scale and the upper bound of the domain.
+        strike = bound(strike, DESIRED_PRECISION_SCALAR + 1, 1e27 - 2); // bounds between the lowest value on the desired precision scale and the upper bound of the domain.
+        v = bound(v, 1, 1e7); // bps
+        t = bound(t, 1, 500 days); // seconds
+
+        // scales v up to ray since we bound by a small range in units of bps, since this is a realistic range
+        v = v * 1e27 / 1e4; // Divides by BPS units which is 1e4, since 1e4 == 100% == 1.0 == 1e27 ray.
+
+        int256 k = 0; // invariant
+        // Assume we want to use `ndtr()` but have an input with WAD units (1E18).
+        // The input can be scaled to match the units of 1E27.
+        // But do we lose precision in the output?
+        uint256 outputWith1E27Input = RMM01.getY(x1, strike, v, t, k);
+        x1 = x1 / 1e9; // Scale down to 1E18.
+        x1 = x1 * 1e9; // Scale back up to 1E27, losing all precision past the 1E18 radix point.
+        uint256 outputWith1E18Input = RMM01.getY(x1, strike, v, t, k);
+
+        // Asserts that the outputs scaled to 1E18 are equal up to 1E18 precision, regardless of the input precision.
+        assertApproxEqPrecision(outputWith1E27Input, outputWith1E18Input, 1e18, "getY-input-precision");
+    }
 }
diff --git a/src/test/RayMathLib.t.sol b/src/test/RayMathLib.t.sol
index 0e78b4c..684ef3c 100644
--- a/src/test/RayMathLib.t.sol
+++ b/src/test/RayMathLib.t.sol
@@ -6,6 +6,7 @@ import "solmate/utils/FixedPointMathLib.sol";
 import "../libraries/RayMathLib.sol" as RayMathLib;
 
 contract TestRayMathLib is Test {
+    /* // todo: once we have a better sqrt.
     function test_logfp() public {
         int256 x = 3.5e27;
         int256 y = RayMathLib.logfp(x);
@@ -13,6 +14,7 @@ contract TestRayMathLib is Test {
         assertEq(y, 1.25276296849536799568812062e27, "logfp-3.5e27");
     }
 
+    // todo: once we have a better sqrt.
     function testFuzz_sqrtfp(int256 x) public {
         vm.assume(x > 0);
 
@@ -20,5 +22,5 @@ contract TestRayMathLib is Test {
         int256 y2 = int256(FixedPointMathLib.sqrt(uint256(x))); // Units of 1e9
 
         assertEq(y / 1e14, y2, "sqrtfp");
-    }
+    } */
 }

From a637c25b278f53eaf55226e97a7e99452ff977a7 Mon Sep 17 00:00:00 2001
From: alex <alexangel.j@gmail.com>
Date: Mon, 22 May 2023 18:31:00 -0700
Subject: [PATCH 9/9] test(extended): moves extended assertions into dedicated
 util

---
 src/test/Ndtr.t.sol                       |  39 +-----
 src/test/Ndtri.t.sol                      |  39 +-----
 src/test/RMM01Monotonicity.t.sol          | 102 +-------------
 src/test/utils/ExtendedAssertionsTest.sol | 154 ++++++++++++++++++++++
 4 files changed, 160 insertions(+), 174 deletions(-)
 create mode 100644 src/test/utils/ExtendedAssertionsTest.sol

diff --git a/src/test/Ndtr.t.sol b/src/test/Ndtr.t.sol
index 4ec48ec..bf22046 100644
--- a/src/test/Ndtr.t.sol
+++ b/src/test/Ndtr.t.sol
@@ -1,47 +1,12 @@
 // SPDX-License-Identifier: GPL-3.0-only
 pragma solidity ^0.8.4;
 
-import "forge-std/Test.sol";
+import "../test/utils/ExtendedAssertionsTest.sol";
 import "../Ndtr.sol";
 import "../Gaussian.sol";
 
 /// @dev Hardcoded expected values computed manually using `normalcdlower` operation in https://keisan.casio.com/calculator
-contract TestNdtr is Test {
-    uint256 constant MAX_ROUNDING_DELTA = 1;
-    int256 constant NDTR_ACCURACY = 1e18;
-
-    /// @notice Compares two in256 values up to a precision with a base of RAY.
-    /// @dev IMPORTANT! Asserts `a` and `b` are within 1 wei up to `precision`.
-    function assertApproxEqPrecision(int256 a, int256 b, int256 precision, string memory message) internal {
-        // Amount to divide by to scale to precision.
-        int256 scalar = RAY / precision;
-
-        // Gets the digits passed the precision end point.
-        uint256 remainder0 = mulmod(uint256(a), uint256(precision), uint256(RAY));
-        uint256 remainder1 = mulmod(uint256(b), uint256(precision), uint256(RAY));
-
-        // For debugging...
-        if (false) {
-            console.log("===== RAW AMOUNTS =====");
-            console.logInt(a);
-            console.logInt(b);
-
-            console.log("===== SCALED AMOUNTS =====");
-            console.logInt(a / scalar);
-            console.logInt(b / scalar);
-
-            console.log("===== REMAINDERS =====");
-            console.log("remainder0", remainder0);
-            console.log("remainder1", remainder1);
-        }
-
-        // Converts units to precision.
-        a = a * precision / RAY;
-        b = b * precision / RAY;
-
-        assertApproxEqAbs(a, b, MAX_ROUNDING_DELTA, message);
-    }
-
+contract TestNdtr is ExtendedAssertionsTest {
     function test_ndtr() public {
         assertApproxEqPrecision(Ndtr.ndtr(RAY), int256(0.841344746068542948585232545e27), NDTR_ACCURACY, "ndtr-1");
         assertApproxEqPrecision(
diff --git a/src/test/Ndtri.t.sol b/src/test/Ndtri.t.sol
index 89ad6c0..017401c 100644
--- a/src/test/Ndtri.t.sol
+++ b/src/test/Ndtri.t.sol
@@ -1,47 +1,12 @@
 // SPDX-License-Identifier: GPL-3.0-only
 pragma solidity ^0.8.4;
 
-import "forge-std/Test.sol";
+import "../test/utils/ExtendedAssertionsTest.sol";
 import "../Ndtri.sol";
 import "../Gaussian.sol";
 
 /// @dev Hardcoded expected values computed manually using `normalcdlower` operation in https://keisan.casio.com/calculator
-contract TestNdtri is Test {
-    int256 constant NDTRI_ACCURACY = 1e18;
-    uint256 constant MAX_ROUNDING_DELTA = 1;
-
-    /// @notice Compares two in256 values up to a precision with a base of RAY.
-    /// @dev IMPORTANT! Asserts `a` and `b` are within 1 wei up to `precision`.
-    function assertApproxEqPrecision(int256 a, int256 b, int256 precision, string memory message) internal {
-        // Amount to divide by to scale to precision.
-        int256 scalar = RAY / precision;
-
-        // Gets the digits passed the precision end point.
-        uint256 remainder0 = mulmod(uint256(a), uint256(precision), uint256(RAY));
-        uint256 remainder1 = mulmod(uint256(b), uint256(precision), uint256(RAY));
-
-        // For debugging...
-        if (false) {
-            console.log("===== RAW AMOUNTS =====");
-            console.logInt(a);
-            console.logInt(b);
-
-            console.log("===== SCALED AMOUNTS =====");
-            console.logInt(a / scalar);
-            console.logInt(b / scalar);
-
-            console.log("===== REMAINDERS =====");
-            console.log("remainder0", remainder0);
-            console.log("remainder1", remainder1);
-        }
-
-        // Converts units to precision.
-        a = a * precision / RAY;
-        b = b * precision / RAY;
-
-        assertApproxEqAbs(a, b, MAX_ROUNDING_DELTA, message);
-    }
-
+contract TestNdtri is ExtendedAssertionsTest {
     /// @dev Tests various inputs on the domain (0, 1), where input is a probability between 0 and 1.
     function test_ndtri() public {
         assertApproxEqPrecision(Ndtri.ndtri(RAY / 2), int256(0), NDTRI_ACCURACY, "ndtri-0.5");
diff --git a/src/test/RMM01Monotonicity.t.sol b/src/test/RMM01Monotonicity.t.sol
index 4709f05..628893d 100644
--- a/src/test/RMM01Monotonicity.t.sol
+++ b/src/test/RMM01Monotonicity.t.sol
@@ -1,7 +1,7 @@
 // SPDX-License-Identifier: GPL-3.0-only
 pragma solidity ^0.8.4;
 
-import "forge-std/Test.sol";
+import "../test/utils/ExtendedAssertionsTest.sol";
 import "../Ndtr.sol";
 import "../Ndtri.sol";
 
@@ -57,105 +57,7 @@ library RMM01 {
     }
 }
 
-contract TestRMM01Monotonicity is Test {
-    bool DEBUG = false;
-
-    uint256 constant MAX_ROUNDING_DELTA = 1;
-    uint256 constant DESIRED_PRECISION = 1e18;
-    uint256 constant DESIRED_PRECISION_SCALAR = 1e9; // Desired precision is 1e18, therefore scalar is 1e27 / 1e18 = 1e9
-
-    /// @notice Compares two in256 values up to a precision with a base of RAY.
-    /// @dev IMPORTANT! Asserts `a` and `b` are within 1 wei up to `precision`.
-    function assertApproxEqPrecision(uint256 a, uint256 b, uint256 precision, string memory message) internal {
-        // Amount to divide by to scale to precision.
-        uint256 scalar = uint256(RAY) / precision;
-
-        // Gets the digits passed the precision end point.
-        uint256 remainder0 = mulmod(uint256(a), uint256(precision), uint256(RAY));
-        uint256 remainder1 = mulmod(uint256(b), uint256(precision), uint256(RAY));
-
-        // For debugging...
-        if (false) {
-            console.log("===== RAW AMOUNTS =====");
-            console.log("a", a);
-            console.log("b", b);
-
-            console.log("===== SCALED AMOUNTS =====");
-            console.log("a / scalar", a / scalar);
-            console.log("b / scalar", b / scalar);
-
-            console.log("===== REMAINDERS =====");
-            console.log("remainder0", remainder0);
-            console.log("remainder1", remainder1);
-        }
-
-        // Converts units to precision.
-        a = a * precision / uint256(RAY);
-        b = b * precision / uint256(RAY);
-
-        assertApproxEqAbs(a, b, MAX_ROUNDING_DELTA, message);
-    }
-
-    /// @dev Asserts `a` is greater than or equal to `b` up to `precision` decimal places.
-    /// @param a The expected larger value.
-    /// @param b The expected smaller value.
-    /// @param precision The number of decimal places to check.
-    /// @param message The message to display if the assertion fails.
-    function assertGtePrecisionNotStrict(uint256 a, uint256 b, uint256 precision, string memory message) internal {
-        // Amount to divide by to scale to precision.
-        uint256 scalar = uint256(RAY) / precision;
-
-        // Quantities beyond the radix point at `precision`. These values are truncated before being checked in this assertion.
-        uint256 remainder0 = mulmod(uint256(a), uint256(precision), uint256(RAY));
-        uint256 remainder1 = mulmod(uint256(b), uint256(precision), uint256(RAY));
-
-        if (DEBUG) {
-            console.log("===== RAW AMOUNTS =====");
-            console.log("a", a);
-            console.log("b", b);
-
-            console.log("===== SCALED AMOUNTS =====");
-            console.log("a", a / scalar);
-            console.log("b", b / scalar);
-
-            console.log("===== REMAINDERS =====");
-            console.log("remainder0", remainder0);
-            console.log("remainder1", remainder1);
-        }
-
-        assertTrue(a / scalar >= b / scalar, "a >= b");
-    }
-
-    /// @dev Asserts `b` is greater than or equal to `a` up to `precision` decimal places.
-    /// @param a The expected smaller value.
-    /// @param b The expected larger value.
-    /// @param precision The number of decimal places to check.
-    /// @param message The message to display if the assertion fails.
-    function assertLtePrecisionNotStrict(uint256 a, uint256 b, uint256 precision, string memory message) internal {
-        // Amount to divide by to scale to precision.
-        uint256 scalar = uint256(RAY) / precision;
-
-        // Quantities beyond the radix point at `precision`. These values are truncated before being checked in this assertion.
-        uint256 remainder0 = mulmod(uint256(a), uint256(precision), uint256(RAY));
-        uint256 remainder1 = mulmod(uint256(b), uint256(precision), uint256(RAY));
-
-        if (DEBUG) {
-            console.log("===== RAW AMOUNTS =====");
-            console.log("a", a);
-            console.log("b", b);
-
-            console.log("===== SCALED AMOUNTS =====");
-            console.log("a", a / scalar);
-            console.log("b", b / scalar);
-
-            console.log("===== REMAINDERS =====");
-            console.log("remainder0", remainder0);
-            console.log("remainder1", remainder1);
-        }
-
-        assertTrue(a / scalar <= b / scalar, "a <= b");
-    }
-
+contract TestRMM01Monotonicity is ExtendedAssertionsTest {
     /// ("✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨");
     /// ("✨✨INVARIANT HAS DECREASING MONOTONICITY✨✨");
     /// ("✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨");
diff --git a/src/test/utils/ExtendedAssertionsTest.sol b/src/test/utils/ExtendedAssertionsTest.sol
new file mode 100644
index 0000000..21c11dd
--- /dev/null
+++ b/src/test/utils/ExtendedAssertionsTest.sol
@@ -0,0 +1,154 @@
+// SPDX-License-Identifier: GPL-3.0-only
+pragma solidity ^0.8.4;
+
+import "forge-std/Test.sol";
+
+contract ExtendedAssertionsTest is Test {
+    ////////////////
+    // Scale      //
+    ////////////////
+    int256 constant SCALE = 1e27;
+    uint256 constant DESIRED_PRECISION = 1e18;
+    uint256 constant DESIRED_PRECISION_SCALAR = 1e9; // Desired precision is 1e18, therefore scalar is 1e27 / 1e18 = 1e9
+
+    ////////////////
+    // Accuracy   //
+    ////////////////
+    int256 constant NDTR_ACCURACY = 1e18;
+    int256 constant NDTRI_ACCURACY = 1e18;
+
+    ////////////////
+    // ApproxEq   //
+    ////////////////
+    uint256 constant MAX_ROUNDING_DELTA = 1;
+
+    ////////////////
+    // Functions  //
+    ////////////////
+
+    /// @notice Compares two in256 values up to a precision with a base of SCALE.
+    /// @dev IMPORTANT! Asserts `a` and `b` are within 1 wei up to `precision`.
+    function assertApproxEqPrecision(uint256 a, uint256 b, uint256 precision, string memory message) internal {
+        // Amount to divide by to scale to precision.
+        uint256 scalar = uint256(SCALE) / precision;
+
+        // Gets the digits passed the precision end point.
+        uint256 remainder0 = mulmod(uint256(a), uint256(precision), uint256(SCALE));
+        uint256 remainder1 = mulmod(uint256(b), uint256(precision), uint256(SCALE));
+
+        // For debugging...
+        if (false) {
+            console.log("===== RAW AMOUNTS =====");
+            console.log("a", a);
+            console.log("b", b);
+
+            console.log("===== SCALED AMOUNTS =====");
+            console.log("a / scalar", a / scalar);
+            console.log("b / scalar", b / scalar);
+
+            console.log("===== REMAINDERS =====");
+            console.log("remainder0", remainder0);
+            console.log("remainder1", remainder1);
+        }
+
+        // Converts units to precision.
+        a = a * precision / uint256(SCALE);
+        b = b * precision / uint256(SCALE);
+
+        assertApproxEqAbs(a, b, MAX_ROUNDING_DELTA, message);
+    }
+
+    /// @notice Compares two in256 values up to a precision with a base of SCALE.
+    /// @dev IMPORTANT! Asserts `a` and `b` are within 1 wei up to `precision`.
+    function assertApproxEqPrecision(int256 a, int256 b, int256 precision, string memory message) internal {
+        // Amount to divide by to scale to precision.
+        int256 scalar = SCALE / precision;
+
+        // Gets the digits passed the precision end point.
+        uint256 remainder0 = mulmod(uint256(a), uint256(precision), uint256(SCALE));
+        uint256 remainder1 = mulmod(uint256(b), uint256(precision), uint256(SCALE));
+
+        // For debugging...
+        if (false) {
+            console.log("===== RAW AMOUNTS =====");
+            console.logInt(a);
+            console.logInt(b);
+
+            console.log("===== SCALED AMOUNTS =====");
+            console.logInt(a / scalar);
+            console.logInt(b / scalar);
+
+            console.log("===== REMAINDERS =====");
+            console.log("remainder0", remainder0);
+            console.log("remainder1", remainder1);
+        }
+
+        // Converts units to precision.
+        a = a * precision / SCALE;
+        b = b * precision / SCALE;
+
+        assertApproxEqAbs(a, b, MAX_ROUNDING_DELTA, message);
+    }
+
+    /// @dev Asserts `a` is greater than or equal to `b` up to `precision` decimal places.
+    /// @param a The expected larger value.
+    /// @param b The expected smaller value.
+    /// @param precision The number of decimal places to check.
+    /// @param message The message to display if the assertion fails.
+    function assertGtePrecisionNotStrict(uint256 a, uint256 b, uint256 precision, string memory message) internal {
+        // Amount to divide by to scale to precision.
+        uint256 scalar = uint256(SCALE) / precision;
+
+        // Quantities beyond the radix point at `precision`. These values are truncated before being checked in this assertion.
+        uint256 remainder0 = mulmod(uint256(a), uint256(precision), uint256(SCALE));
+        uint256 remainder1 = mulmod(uint256(b), uint256(precision), uint256(SCALE));
+
+        // For debugging...
+        if (false) {
+            console.log("===== RAW AMOUNTS =====");
+            console.log("a", a);
+            console.log("b", b);
+
+            console.log("===== SCALED AMOUNTS =====");
+            console.log("a", a / scalar);
+            console.log("b", b / scalar);
+
+            console.log("===== REMAINDERS =====");
+            console.log("remainder0", remainder0);
+            console.log("remainder1", remainder1);
+        }
+
+        assertTrue(a / scalar >= b / scalar, "a >= b");
+    }
+
+    /// @dev Asserts `b` is greater than or equal to `a` up to `precision` decimal places.
+    /// @param a The expected smaller value.
+    /// @param b The expected larger value.
+    /// @param precision The number of decimal places to check.
+    /// @param message The message to display if the assertion fails.
+    function assertLtePrecisionNotStrict(uint256 a, uint256 b, uint256 precision, string memory message) internal {
+        // Amount to divide by to scale to precision.
+        uint256 scalar = uint256(SCALE) / precision;
+
+        // Quantities beyond the radix point at `precision`. These values are truncated before being checked in this assertion.
+        uint256 remainder0 = mulmod(uint256(a), uint256(precision), uint256(SCALE));
+        uint256 remainder1 = mulmod(uint256(b), uint256(precision), uint256(SCALE));
+
+        // For debugging...
+        if (false) {
+            console.log("===== RAW AMOUNTS =====");
+            console.log("a", a);
+            console.log("b", b);
+
+            console.log("===== SCALED AMOUNTS =====");
+            console.log("a", a / scalar);
+            console.log("b", b / scalar);
+
+            console.log("===== REMAINDERS =====");
+            console.log("remainder0", remainder0);
+            console.log("remainder1", remainder1);
+        }
+
+        assertTrue(a / scalar <= b / scalar, "a <= b");
+    }
+}