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
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/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/Ndtr.sol b/src/Ndtr.sol
new file mode 100644
index 0000000..8180a7f
--- /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 = expfp(z);
+
+        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 = x * 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..191e81f
--- /dev/null
+++ b/src/Ndtri.sol
@@ -0,0 +1,198 @@
+// 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();
+
+    ////////////////
+    // Constants  //
+    ////////////////
+
+    int256 constant RAY_EXPNEG2 = 0.135335283236612691893999494e27;
+
+    //////////////////////////////////////////////
+    // 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;
+
+        // Use a different approximation if y > exp(-2).
+        if (y > RAY_ONE - RAY_EXPNEG2) {
+            // 0.135... = exp(-2)
+            y = RAY_ONE - y;
+            code = 0;
+        }
+
+        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;
+            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.
+
+        // 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)
+        // 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 - x1;
+        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..1ff86b0
--- /dev/null
+++ b/src/libraries/RayMathLib.sol
@@ -0,0 +1,78 @@
+// 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) 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) pure returns (int256) {
+    x = x / 1e9; // Scale from 1e27 -> 1e18.
+
+    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
+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
new file mode 100644
index 0000000..bf22046
--- /dev/null
+++ b/src/test/Ndtr.t.sol
@@ -0,0 +1,266 @@
+// SPDX-License-Identifier: GPL-3.0-only
+pragma solidity ^0.8.4;
+
+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 ExtendedAssertionsTest {
+    function test_ndtr() public {
+        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"
+        );
+        assertApproxEqPrecision(
+            Ndtr.ndtr(RAY / 100), int256(0.503989356314631603924543868e27), NDTR_ACCURACY, "ndtr-0.01"
+        );
+        // negative
+        assertApproxEqPrecision(
+            Ndtr.ndtr(-RAY / 10000000), int256(0.499999960105771959856203402e27), NDTR_ACCURACY, "ndtr--0.0000001"
+        );
+        assertApproxEqPrecision(
+            Ndtr.ndtr(-RAY / 100), int256(0.496010643685368396075456132e27), NDTR_ACCURACY, "ndtr--0.01"
+        );
+        // 0
+        assertApproxEqPrecision(Ndtr.ndtr(0), int256(0.5e27), NDTR_ACCURACY, "ndtr-0");
+        // random values at scale of 1e27
+        assertApproxEqPrecision(
+            Ndtr.ndtr(0.123456789e27), int256(0.5491273050781420888711e27), NDTR_ACCURACY, "ndtr-0.123456789"
+        );
+        assertApproxEqPrecision(
+            Ndtr.ndtr(0.987654321e27), int256(0.83833901356624443490786e27), NDTR_ACCURACY, "ndtr-0.987654321"
+        );
+        assertApproxEqPrecision(
+            Ndtr.ndtr(-0.123456789e27), int256(0.4508726949218579111288e27), NDTR_ACCURACY, "ndtr--0.123456789"
+        );
+        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 {
+        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"
+        );
+        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"
+        );
+        assertApproxEqPrecision(
+            Ndtr.ndtr(RAY / NDTR_ACCURACY * 1e1), int256(0.500000000000000003989422804e27), NDTR_ACCURACY, "ndtr-1e-17"
+        );
+        assertApproxEqPrecision(
+            Ndtr.ndtr(-RAY / NDTR_ACCURACY * 1e1),
+            int256(0.499999999999999996010577195e27),
+            NDTR_ACCURACY,
+            "ndtr--1e-17"
+        );
+        assertApproxEqPrecision(
+            Ndtr.ndtr(RAY / NDTR_ACCURACY * 1e2), int256(0.50000000000000003989422804e27), NDTR_ACCURACY, "ndtr-1e-16"
+        );
+        assertApproxEqPrecision(
+            Ndtr.ndtr(-RAY / NDTR_ACCURACY * 1e2),
+            int256(0.499999999999999960105771959e27),
+            NDTR_ACCURACY,
+            "ndtr--1e-16"
+        );
+        assertApproxEqPrecision(
+            Ndtr.ndtr(RAY / NDTR_ACCURACY * 1e3), int256(0.500000000000000398942280401e27), NDTR_ACCURACY, "ndtr-1e-15"
+        );
+        assertApproxEqPrecision(
+            Ndtr.ndtr(-RAY / NDTR_ACCURACY * 1e3),
+            int256(0.499999999999999601057719598e27),
+            NDTR_ACCURACY,
+            "ndtr--1e-15"
+        );
+        assertApproxEqPrecision(
+            Ndtr.ndtr(RAY / NDTR_ACCURACY * 1e4), int256(0.500000000000003989422804014e27), NDTR_ACCURACY, "ndtr-1e-14"
+        );
+        assertApproxEqPrecision(
+            Ndtr.ndtr(-RAY / NDTR_ACCURACY * 1e4),
+            int256(0.499999999999996010577195985e27),
+            NDTR_ACCURACY,
+            "ndtr--1e-14"
+        );
+        assertApproxEqPrecision(
+            Ndtr.ndtr(RAY / NDTR_ACCURACY * 1e5), int256(0.500000000000039894228040143e27), NDTR_ACCURACY, "ndtr-1e-13"
+        );
+        assertApproxEqPrecision(
+            Ndtr.ndtr(-RAY / NDTR_ACCURACY * 1e5),
+            int256(0.499999999999960105771959856e27),
+            NDTR_ACCURACY,
+            "ndtr--1e-13"
+        );
+    }
+
+    function test_ndtr_far_zero() public {
+        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
+        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 {
+        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"
+        );
+        assertApproxEqPrecision(
+            Ndtr.ndtr(RAY * 2 / 3), int256(0.747507462453077086935938175e27), NDTR_ACCURACY, "ndtr-0.666666"
+        );
+        assertApproxEqPrecision(
+            Ndtr.ndtr(RAY * 5 / 6), int256(0.7976716190363569746314664e27), NDTR_ACCURACY, "ndtr-0.833333"
+        );
+        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"
+        );
+        assertApproxEqPrecision(
+            Ndtr.ndtr(-RAY * 2 / 3), int256(0.252492537546922913064061824e27), NDTR_ACCURACY, "ndtr--0.666666"
+        );
+        assertApproxEqPrecision(
+            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
+    /// @dev As x -> inf, y -> 1
+    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
+    /// @dev As x -> -inf, y -> 0
+    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?
+    }
+
+    /// 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);
+        vm.assume(z < -RAY_MAXLOG);
+        int256 y = Ndtr.erfc(x);
+        assertEq(y, int256(RAY_TWO), "erfc-not-two");
+    }
+
+    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");
+    }
+
+    /// @dev erfc domain is -2 < 0 < 2, check between (-2, 0)
+    function testFuzz_erfc_NegativeInputIsBounded(int256 x) public {
+        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 > 1e9);
+        vm.assume(x < RAY_TWO);
+        int256 y = Ndtr.erfc(x);
+        assertGe(y, 0 ether);
+        assertLe(y, RAY);
+    }
+
+    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/Ndtri.t.sol b/src/test/Ndtri.t.sol
new file mode 100644
index 0000000..017401c
--- /dev/null
+++ b/src/test/Ndtri.t.sol
@@ -0,0 +1,160 @@
+// SPDX-License-Identifier: GPL-3.0-only
+pragma solidity ^0.8.4;
+
+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 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");
+        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"
+        );
+        assertApproxEqPrecision(
+            Ndtri.ndtri(RAY * 77 / 93), int256(0.946122671364591906727098902e27), NDTRI_ACCURACY, "ndtri-0.827956"
+        );
+        assertApproxEqPrecision(
+            Ndtri.ndtri(RAY * 11 / 85), int256(-1.1291761577077979287868872e27), NDTRI_ACCURACY, "ndtri-0.129412"
+        );
+
+        assertApproxEqPrecision(
+            Ndtri.ndtri(729329433526774616220000000),
+            int256(0.610786196335546300516717471e27),
+            NDTRI_ACCURACY,
+            "ndtri-0.729329"
+        );
+    }
+
+    /// @dev Starting at 1e27, decrements by atleast RAY / NDTRI_ACCURACY = 1e9, since thats the lowest unit
+    /// that has precision.
+    /// 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;
+        assertApproxEqPrecision(
+            Ndtri.ndtri(RAY - decrement),
+            int256(8.75729034878231506388112862e27),
+            NDTRI_ACCURACY,
+            "ndtri-0.999999999999999999000000000"
+        );
+        assertApproxEqPrecision(
+            Ndtri.ndtri(RAY - decrement * 1e1), int256(8.49379322410959807444471881e27), NDTRI_ACCURACY, "ndtri-1-1e-17"
+        );
+        assertApproxEqPrecision(
+            Ndtri.ndtri(RAY - decrement * 1e2), int256(8.22208221613043561267585878e27), NDTRI_ACCURACY, "ndtri-1-1e-16"
+        );
+        assertApproxEqPrecision(
+            Ndtri.ndtri(RAY - decrement * 1e3), int256(7.94134532617099678096674357e27), NDTRI_ACCURACY, "ndtri-1-1e-15"
+        );
+        assertApproxEqPrecision(
+            Ndtri.ndtri(RAY - decrement * 1e4), int256(7.65062809293526881641896581e27), NDTRI_ACCURACY, "ndtri-1-1e-14"
+        );
+        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
+    /// 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"
+        );
+        assertApproxEqPrecision(
+            Ndtri.ndtri(RAY / NDTRI_ACCURACY * 1e1), -int256(8.493793224109598e27), NDTRI_ACCURACY, "ndtri-1e-17"
+        );
+        assertApproxEqPrecision(
+            Ndtri.ndtri(RAY / NDTRI_ACCURACY * 1e2), -int256(8.2220822161304356e27), NDTRI_ACCURACY, "ndtri-1e-16"
+        );
+        assertApproxEqPrecision(
+            Ndtri.ndtri(RAY / NDTRI_ACCURACY * 1e3), -int256(7.9413453261709968e27), NDTRI_ACCURACY, "ndtri-1e-15"
+        );
+        assertApproxEqPrecision(
+            Ndtri.ndtri(RAY / NDTRI_ACCURACY * 1e4), -int256(7.65062809293526882e27), NDTRI_ACCURACY, "ndtri-1e-14"
+        );
+        assertApproxEqPrecision(
+            Ndtri.ndtri(RAY / NDTRI_ACCURACY * 1e5), -int256(7.348796102800677518e27), NDTRI_ACCURACY, "ndtri-1e-13"
+        );
+    }
+
+    /// @dev Expected `b` values computed manually using `normalcdlower` operation in https://keisan.casio.com/calculator
+    /// 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..."
+        );
+        assertApproxEqPrecision(
+            Ndtri.ndtri(RAY * 3 / 5), int256(0.253347103135799798798196181e27), NDTRI_ACCURACY, "ndtri-0.6"
+        );
+        assertApproxEqPrecision(
+            Ndtri.ndtri(RAY * 2 / 3), int256(0.43072729929545749020594039e27), NDTRI_ACCURACY, "ndtri-0.666666..."
+        );
+        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: 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);
+
+        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!
+    /// @dev As x -> 0, ndtri(x) -> -inf
+    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");
+    }
+
+    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/RMM01Monotonicity.t.sol b/src/test/RMM01Monotonicity.t.sol
new file mode 100644
index 0000000..628893d
--- /dev/null
+++ b/src/test/RMM01Monotonicity.t.sol
@@ -0,0 +1,138 @@
+// SPDX-License-Identifier: GPL-3.0-only
+pragma solidity ^0.8.4;
+
+import "../test/utils/ExtendedAssertionsTest.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;
+
+    /// @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
+            int256 sqrtTauVol = mulfp(int256(volatilityRay), sqrtTau);
+            int256 phi = int256(PRECISION) - int256(x);
+            int256 inverseCdf = Ndtri.ndtri(phi);
+            int256 input = inverseCdf - sqrtTauVol;
+            int256 cdf = Ndtr.ndtr(input);
+            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 ExtendedAssertionsTest {
+    /// ("✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨");
+    /// ("✨✨INVARIANT HAS DECREASING MONOTONICITY✨✨");
+    /// ("✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨");
+    /// As x -> 1, y -> 0
+    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
+        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
+        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);
+        console.log("y1", y1);
+        console.log("x2", x1 + 1);
+        console.log("y2", y2);
+
+        // 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 INCREASING MONOTONICITY✨✨");
+    /// ("✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨✨");
+    /// As x -> 0, y -> K
+    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
+        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
+        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);
+        console.log("y1", y1);
+        console.log("x2", x1 - 1);
+        console.log("y2", y2);
+
+        // As x gets smaller, expect y to get larger.
+        // 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
new file mode 100644
index 0000000..684ef3c
--- /dev/null
+++ b/src/test/RayMathLib.t.sol
@@ -0,0 +1,26 @@
+// 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 {
+    /* // todo: once we have a better sqrt.
+    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");
+    }
+
+    // todo: once we have a better sqrt.
+    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");
+    } */
+}
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");
+    }
+}