diff --git a/examples/toffoli.cpp b/examples/toffoli.cpp
new file mode 100644
index 000000000..267e5feaa
--- /dev/null
+++ b/examples/toffoli.cpp
@@ -0,0 +1,76 @@
+// Toffoli gate simulation
+// Source: ./examples/toffoli.cpp
+
+#include <iostream>
+
+#include "qpp/qpp.h"
+
+int main() {
+    using namespace qpp;
+
+    std::cout << ">> Toffoli gate simulation\n\n";
+
+    ket psi_in = randket(8);
+    std::cout << ">> Input state:\n";
+    std::cout << disp(dirac(psi_in)) << "\n\n";
+
+    /**
+     * Toffoli gate (control control not)
+     *
+     *  ---+---
+     *     |
+     *  ---+---
+     *     |
+     *  ---X---
+     */
+    ket result = apply(psi_in, gt.TOF, {0, 1, 2});
+    std::cout << ">> Toffoli gate output state:\n";
+    std::cout << disp(dirac(result)) << "\n\n";
+
+    /**
+     * Toffoli with T and CNOT
+     *
+     *  -------------+-------------+-----+---T---+--
+     *               |             |     |       |
+     *  -----+-------------+----------T--X--T_d--X--
+     *       |       |     |       |
+     *  --H--X--T_d--X--T--X--T_d--X--T--H----------
+     */
+    result = apply(psi_in, gt.H, {2});
+    result = applyCTRL(result, gt.X, {1}, {2});
+    result = apply(result, adjoint(gt.T), {2});
+    result = applyCTRL(result, gt.X, {0}, {2});
+    result = apply(result, gt.T, {2});
+    result = applyCTRL(result, gt.X, {1}, {2});
+    result = apply(result, adjoint(gt.T), {2});
+    result = applyCTRL(result, gt.X, {0}, {2});
+    result = apply(result, gt.T, {1});
+    result = apply(result, gt.T, {2});
+    result = applyCTRL(result, gt.X, {0}, {1});
+    result = apply(result, gt.T, {0});
+    result = apply(result, adjoint(gt.T), {1});
+    result = apply(result, gt.H, {2});
+    result = applyCTRL(result, gt.X, {0}, {1});
+    std::cout << ">> Toffoli with T and CNOT output state:\n";
+    std::cout << disp(dirac(result)) << "\n\n";
+
+    /**
+     * Sleator Weinfurter construction
+     * V * V = X
+     *
+     *  -----+-------+---+---
+     *       |       |   |
+     *  --+--X---+---X-------
+     *    |      |       |
+     *  --V-----V_d------V---
+     */
+    cmat sqrtx{cmat::Zero(2, 2)};
+    sqrtx << 0.5 + 0.5 * 1_i, 0.5 - 0.5 * 1_i, 0.5 - 0.5 * 1_i, 0.5 + 0.5 * 1_i;
+    result = applyCTRL(psi_in, sqrtx, {1}, {2});
+    result = applyCTRL(result, gt.X, {0}, {1});
+    result = applyCTRL(result, adjoint(sqrtx), {1}, {2});
+    result = applyCTRL(result, gt.X, {0}, {1});
+    result = applyCTRL(result, sqrtx, {0}, {2});
+    std::cout << ">> Sleator Weinfurter construction output state:\n";
+    std::cout << disp(dirac(result)) << "\n\n";
+}