diff --git a/gaussian b/gaussian
new file mode 100755
index 0000000..7694ef6
Binary files /dev/null and b/gaussian differ
diff --git a/src/Gaussian.sol b/src/Gaussian.sol
index b31d5af..bfbc534 100644
--- a/src/Gaussian.sol
+++ b/src/Gaussian.sol
@@ -27,6 +27,7 @@ library Gaussian {
     error Infinity();
     error NegativeInfinity();
     error Overflow();
+    error OutOfBounds();
 
     uint256 internal constant HALF_WAD = 0.5 ether;
     uint256 internal constant PI = 3_141592653589793238;
@@ -77,6 +78,10 @@ library Gaussian {
      * @custom:source https://mathworld.wolfram.com/Erfc.html.
      */
     function erfc(int256 input) internal pure returns (int256 output) {
+        if (input == 0) {
+            return 1 ether;
+        }
+
         uint256 z = abs(input);
         int256 t;
         int256 step;
@@ -171,6 +176,11 @@ library Gaussian {
      * @custom:source https://mathworld.wolfram.com/InverseErfc.html.
      */
     function ierfc(int256 x) internal pure returns (int256 z) {
+        if (x == 0 || x == 2 ether) revert Infinity();
+        if (x < 0 || x > 2 ether) revert OutOfBounds();
+
+        if (x == SCALAR) return 0;
+
         assembly {
             // x >= 2, iszero(x < 2 ? 1 : 0) ? 1 : 0.
             if iszero(slt(x, TWO)) {
@@ -207,7 +217,14 @@ library Gaussian {
         int256 r;
         assembly {
             function muli(pxn, pxd) -> res {
-                res := sdiv(mul(pxn, pxd), ONE)
+                res := mul(pxn, pxd)
+
+                if iszero(eq(sdiv(res, pxn), pxd)) {
+                    mstore(0, 0x35278d1200000000000000000000000000000000000000000000000000000000)
+                    revert(0, 4)
+                }
+
+                res := sdiv(res, ONE)
             }
 
             r := muli(
@@ -275,12 +292,14 @@ library Gaussian {
      */
     function cdf(int256 x) internal pure returns (int256 z) {
         int256 negated;
+
         assembly {
             let res := sdiv(mul(x, ONE), SQRT2)
             negated := add(not(res), 1)
         }
 
         int256 _erfc = erfc(negated);
+
         assembly {
             z := sdiv(mul(ONE, _erfc), TWO)
         }
@@ -297,9 +316,11 @@ library Gaussian {
      */
     function pdf(int256 x) internal pure returns (int256 z) {
         int256 e;
+
         assembly {
             e := sdiv(mul(add(not(x), 1), x), TWO) // (-x * x) / 2.
         }
+
         e = FixedPointMathLib.expWad(e);
 
         assembly {
diff --git a/src/test/Cdf.t.sol b/src/test/Cdf.t.sol
new file mode 100644
index 0000000..322accf
--- /dev/null
+++ b/src/test/Cdf.t.sol
@@ -0,0 +1,29 @@
+// SPDX-License-Identifier: MIT
+pragma solidity 0.8.13;
+
+import "forge-std/Test.sol";
+
+import {Gaussian} from "../Gaussian.sol";
+
+contract TestCdf is Test {
+    function testDiff_cdf(int256 x) public {
+        vm.assume(x > -2828427124746190093171572875253809907);
+        vm.assume(x < 2828427124746190093171572875253809907);
+        string[] memory inputs = new string[](3);
+        inputs[0] = "./gaussian";
+        inputs[1] = "cdf";
+        inputs[2] = vm.toString(x);
+        bytes memory res = vm.ffi(inputs);
+        uint256 ref = abi.decode(res, (uint256));
+        int256 y = Gaussian.cdf(x);
+
+        // When outputs are very small, we tolerate a larger error
+        if (ref < 1_000_000_000 && y < 1_000_000_000) {
+            // 0.1% of difference
+            assertApproxEqRel(ref, uint256(y), 0.001 ether);
+        } else {
+            // 0.00005% of difference
+            assertApproxEqRel(ref, uint256(y), 0.0000005 ether);
+        }
+    }
+}
diff --git a/src/test/Erfc.t.sol b/src/test/Erfc.t.sol
new file mode 100644
index 0000000..546a6e0
--- /dev/null
+++ b/src/test/Erfc.t.sol
@@ -0,0 +1,59 @@
+// SPDX-License-Identifier: MIT
+pragma solidity 0.8.13;
+
+import "forge-std/Test.sol";
+
+import {Gaussian} from "../Gaussian.sol";
+
+contract TestErfc is Test {
+    function testFuzz_erfc_RevertWhenInputIsTooLow(int256 x) public {
+        vm.assume(x <= -1999999999999999998000000000000000002);
+
+        // TODO: Investigate why the error selector is 0x4d2d75b1
+        // instead of 0x35278d12
+        vm.expectRevert();
+        int256 y = Gaussian.erfc(x);
+        y;
+    }
+
+    function testFuzz_erfc_RevertWhenInputIsTooHigh(int256 x) public {
+        vm.assume(x > 1999999999999999998000000000000000002);
+        vm.expectRevert(Gaussian.Overflow.selector);
+        int256 y = Gaussian.erfc(x);
+        y;
+    }
+
+    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);
+    }
+
+    function testDiff_erfc(int256 x) public {
+        vm.assume(x < 1999999999999999998000000000000000002);
+        vm.assume(x > -1999999999999999998000000000000000002);
+        string[] memory inputs = new string[](3);
+        inputs[0] = "./gaussian";
+        inputs[1] = "erfc";
+        inputs[2] = vm.toString(x);
+        bytes memory res = vm.ffi(inputs);
+        uint256 ref = abi.decode(res, (uint256));
+        int256 y = Gaussian.erfc(x);
+        // Results have a 0.0001% difference
+        assertApproxEqRel(ref, uint256(y), 0.000001 ether);
+    }
+}
diff --git a/src/test/Gaussian.t.sol b/src/test/Gaussian.t.sol
index ffbf394..f3eea57 100644
--- a/src/test/Gaussian.t.sol
+++ b/src/test/Gaussian.t.sol
@@ -157,6 +157,7 @@ contract TestGaussian is Test {
     function testReference_erfc_Equality(int256 input) public {
         vm.assume(input < HIGH);
         vm.assume(input > LOW);
+        vm.assume(input != 0);
         int256 actual = Gaussian.erfc(input);
         int256 expected = Ref.erfc(input);
         assertEq(actual, expected, "erfc-inequality");
@@ -164,9 +165,8 @@ contract TestGaussian is Test {
 
     // todo: investigate, reverts with Infinity() if input is 1 because of rounding down?
     function testReference_ierfc_Equality(int256 input) public {
-        vm.assume(input < HIGH);
-        vm.assume(input > LOW);
-        vm.assume(input != 1);
+        vm.assume(input < 2 ether);
+        vm.assume(input > 1);
         int256 actual = Gaussian.ierfc(input);
         int256 expected = Ref.ierfc(input);
         assertEq(actual, expected, "ierfc-inequality");
@@ -175,9 +175,13 @@ contract TestGaussian is Test {
     function testReference_cdf_Equality(int256 input) public {
         vm.assume(input < HIGH);
         vm.assume(input > LOW);
+        vm.assume(input != 0);
         int256 actual = Gaussian.cdf(input);
         int256 expected = Ref.cdf(input);
-        assertEq(actual, expected, "cdf-inequality");
+        console.logInt(actual);
+        console.logInt(expected);
+        // assertEq(actual, expected, "cdf-inequality");
+        assertApproxEqRel(actual, expected, 0.0015 ether);
     }
 
     function testReference_pdf_Equality(int256 input) public {
diff --git a/src/test/Ierfc.t.sol b/src/test/Ierfc.t.sol
new file mode 100644
index 0000000..cdcae82
--- /dev/null
+++ b/src/test/Ierfc.t.sol
@@ -0,0 +1,54 @@
+// SPDX-License-Identifier: MIT
+pragma solidity 0.8.13;
+
+import "forge-std/Test.sol";
+
+import {Gaussian} from "../Gaussian.sol";
+
+contract TestIerfc is Test {
+    function testFuzz_ierfc_RevertWhenInputIsOutOfBounds(int256 x) public {
+        vm.assume(x < 0 || x > 2 ether);
+        vm.expectRevert(Gaussian.OutOfBounds.selector);
+        int256 y = Gaussian.ierfc(x);
+        y;
+    }
+
+    function test_ierfc_InputOneWillTriggerInfinity() public {
+        vm.expectRevert(Gaussian.Infinity.selector);
+        int256 y = Gaussian.ierfc(1);
+        console.logInt(y);
+    }
+
+    function test_ierfc_ZeroTriggersInfinity() public {
+        vm.expectRevert(Gaussian.Infinity.selector);
+        int256 y = Gaussian.ierfc(0);
+        y;
+    }
+
+    function test_ierfc_TwoTriggersInfinity() public {
+        vm.expectRevert(Gaussian.Infinity.selector);
+        int256 y = Gaussian.ierfc(2 ether);
+        y;
+    }
+
+    function testDiff_ierfc(int64 x) public {
+        vm.assume(x > 0.00001 ether);
+        vm.assume(x < 2 ether);
+        string[] memory inputs = new string[](3);
+        inputs[0] = "./gaussian";
+        inputs[1] = "ierfc";
+        inputs[2] = vm.toString(x);
+        bytes memory res = vm.ffi(inputs);
+        int256 ref = abi.decode(res, (int256));
+        int256 y = Gaussian.ierfc(int256(x));
+
+        // When inputs are very close to 1, we tolerate a larger error
+        if (x > 0.99 ether && x < 1.05 ether) {
+            // 0.15% of difference
+            assertApproxEqRel(ref, y, 0.0015 ether);
+        } else {
+            // 0.0003% of difference
+            assertApproxEqRel(ref, y, 0.000003 ether);
+        }
+    }
+}
diff --git a/src/test/Pdf.t.sol b/src/test/Pdf.t.sol
new file mode 100644
index 0000000..422f1b2
--- /dev/null
+++ b/src/test/Pdf.t.sol
@@ -0,0 +1,29 @@
+// SPDX-License-Identifier: MIT
+pragma solidity 0.8.13;
+
+import "forge-std/Test.sol";
+
+import {Gaussian} from "../Gaussian.sol";
+
+contract TestPdf is Test {
+    function testDiff_pdf(int256 x) public {
+        vm.assume(x > -2828427124746190093171572875253809907);
+        vm.assume(x < 2828427124746190093171572875253809907);
+        string[] memory inputs = new string[](3);
+        inputs[0] = "./gaussian";
+        inputs[1] = "pdf";
+        inputs[2] = vm.toString(x);
+        bytes memory res = vm.ffi(inputs);
+        int256 ref = abi.decode(res, (int256));
+        int256 y = Gaussian.pdf(x);
+
+        // When outputs are very small, we tolerate a larger error
+        if (ref < 1_000_000_000 && y < 1_000_000_000) {
+            // 0.1% of difference
+            assertApproxEqRel(ref, y, 0.001 ether);
+        } else {
+            // 0.00005% of difference
+            assertApproxEqRel(ref, y, 0.0000005 ether);
+        }
+    }
+}
diff --git a/src/test/Ppf.t.sol b/src/test/Ppf.t.sol
new file mode 100644
index 0000000..ae1d596
--- /dev/null
+++ b/src/test/Ppf.t.sol
@@ -0,0 +1,22 @@
+// SPDX-License-Identifier: MIT
+pragma solidity 0.8.13;
+
+import "forge-std/Test.sol";
+
+import {Gaussian} from "../Gaussian.sol";
+
+contract TestPpf is Test {
+    function testDiff_ppf(int64 x) public {
+        vm.assume(x > 0.0000001 ether);
+        vm.assume(x < 1 ether);
+        string[] memory inputs = new string[](3);
+        inputs[0] = "./gaussian";
+        inputs[1] = "ppf";
+        inputs[2] = vm.toString(x);
+        bytes memory res = vm.ffi(inputs);
+        int256 ref = abi.decode(res, (int256));
+        int256 y = Gaussian.ppf(int256(x));
+        // Results have a 0.00165% difference
+        assertApproxEqRel(ref, y, 0.001 ether);
+    }
+}