diff --git a/vm/Cargo.toml b/vm/Cargo.toml
index 84a95c5ffa..7e408c7788 100644
--- a/vm/Cargo.toml
+++ b/vm/Cargo.toml
@@ -15,3 +15,4 @@ serde_derive = "1.0.66"
 serde_json = "1.0.26"
 byteorder = "1.2.6"
 regex = "1"
+statrs = "0.10.0"
diff --git a/vm/src/lib.rs b/vm/src/lib.rs
index f2f61fac73..cfcaf94339 100644
--- a/vm/src/lib.rs
+++ b/vm/src/lib.rs
@@ -15,6 +15,7 @@ extern crate num_complex;
 extern crate num_traits;
 extern crate serde;
 extern crate serde_json;
+extern crate statrs;
 
 extern crate rustpython_parser;
 
diff --git a/vm/src/stdlib/math.rs b/vm/src/stdlib/math.rs
index 0ee716b9bd..0cc800b129 100644
--- a/vm/src/stdlib/math.rs
+++ b/vm/src/stdlib/math.rs
@@ -8,6 +8,8 @@ use super::super::pyobject::{
     DictProtocol, PyContext, PyFuncArgs, PyObjectRef, PyResult, TypeProtocol,
 };
 use super::super::VirtualMachine;
+use statrs::function::erf::{erf, erfc};
+use statrs::function::gamma::{gamma, ln_gamma};
 use std;
 
 // Helper macro:
@@ -156,8 +158,7 @@ fn math_erf(vm: &mut VirtualMachine, args: PyFuncArgs) -> PyResult {
     if x.is_nan() {
         Ok(vm.ctx.new_float(x))
     } else {
-        // TODO: implement algorithm
-        unimplemented!("TODO");
+        Ok(vm.ctx.new_float(erf(x)))
     }
 }
 
@@ -168,8 +169,7 @@ fn math_erfc(vm: &mut VirtualMachine, args: PyFuncArgs) -> PyResult {
     if x.is_nan() {
         Ok(vm.ctx.new_float(x))
     } else {
-        // TODO: implement algorithm
-        unimplemented!("TODO");
+        Ok(vm.ctx.new_float(erfc(x)))
     }
 }
 
@@ -178,8 +178,7 @@ fn math_gamma(vm: &mut VirtualMachine, args: PyFuncArgs) -> PyResult {
     let x = objfloat::get_value(value);
 
     if x.is_finite() {
-        // TODO: implement algorithm
-        unimplemented!("TODO");
+        Ok(vm.ctx.new_float(gamma(x)))
     } else {
         if x.is_nan() || x.is_sign_positive() {
             Ok(vm.ctx.new_float(x))
@@ -194,8 +193,7 @@ fn math_lgamma(vm: &mut VirtualMachine, args: PyFuncArgs) -> PyResult {
     let x = objfloat::get_value(value);
 
     if x.is_finite() {
-        // TODO: implement algorithm
-        unimplemented!("TODO");
+        Ok(vm.ctx.new_float(ln_gamma(x)))
     } else {
         if x.is_nan() {
             Ok(vm.ctx.new_float(x))