diff --git a/include/real/irrational_helpers.hpp b/include/real/irrational_helpers.hpp
index ab49383..248242a 100644
--- a/include/real/irrational_helpers.hpp
+++ b/include/real/irrational_helpers.hpp
@@ -93,49 +93,72 @@ namespace boost {
                 return current_value.digits[n];
             }
 
-            int pi_get_nth_digit(unsigned int n) {
-                if (n > 2000)
-                    throw(pi_precision_exception());
-                // 2000 digits of pi
-                exact_number pi("3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912983367336244065664308602139494639522473719070217986094370277053921717629317675238467481846766940513200056812714526356082778577134275778960917363717872146844090122495343014654958537105079227968925892354201995611212902196086403441815981362977477130996051870721134999999837297804995105973173281609631859502445945534690830264252230825334468503526193118817101000313783875288658753320838142061717766914730359825349042875546873115956286388235378759375195778185778053217122680661300192787661119590921642019893809525720106548586327886593615338182796823030195203530185296899577362259941389124972177528347913151557485724245415069595082953311686172785588907509838175463746493931925506040092770167113900984882401285836160356370766010471018194295559619894676783744944825537977472684710404753464620804668425906949129331367702898915210475216205696602405803815019351125338243003558764024749647326391419927260426992279678235478163600934172164121992458631503028618297455570674983850549458858692699569092721079750930295532116534498720275596023648066549911988183479775356636980742654252786255181841757467289097777279380008164706001614524919217321721477235014144197356854816136115735255213347574184946843852332390739414333454776241686251898356948556209921922218427255025425688767179049460165346680498862723279178608578438382796797668145410095388378636095068006422512520511739298489608412848862694560424196528502221066118630674427862203919494504712371378696095636437191728746776465757396241389086583264599581339047802759009");
-                return pi[n];
+            // TODO optimize: fast binary exponentiation
+            static long mul_exp2_mod(long p, long exp, long m) {
+                long result = p;
+                for(long i = 0; i < exp; i++) {
+                    result = (result * 2) % m;
+                }
+                return result % m;
             }
 
-            // e^x = sum_{k=0}^\inf = x^0/0! + x^1/1! + x^2/2! + x^3/3! + x^3/6! + ...
-
-            class exponential {
-                private:
-                // would be nice to interoperate between long, int, and boost::real::real,
-                // and have ctors from the integral types
-                boost::real::real k_prev = boost::real::real_explicit("0");
-                boost::real::real const * const x_ptr;
-                boost::real::real last_term; // x^kn / kn!
-                boost::real::real current_value; // summation from k0 to k_n, with precision digits
-
-                public:
-                exponential(boost::real::real &x) : x_ptr(&x) {
-                    last_term = boost::real::real ("1");
-                    current_value = boost::real::real ("1");
-                };
-
-                int get_nth_digit(unsigned int n) {
-                    boost::real::real one = boost::real::real_explicit("1");
-                    // if n < k_prev, reset
-
-                    boost::real::real min_bound;
-                    std::get<boost::real::real_explicit>(min_bound.get_real_number()).digits = {1};
-                    std::get<boost::real::real_explicit>(min_bound.get_real_number()).exponent = 1-n;
-
-                    // keep getting terms from the taylor series until the terms go below our precision bound
-                    while((last_term > min_bound) || (last_term == min_bound)) {
-                        last_term = last_term * (*x_ptr) / (k_prev + one);
-                        current_value = current_value + last_term;
-                    }
+            /**
+             * @brief The function returns the n-th digit of binary expansion of pi.
+             * @param n - The number digit index.
+             * @return The n-th digit of pi's binary expansion. 
+             */
 
-                    return 0;
+            int pi_binary_get_nth_digit(unsigned int n) {
+                if (n == 1) return 1;
+                else if (n == 2) return 1;
+                n -= 3; // BBP algorithm computes the (n+1)-th *fractional bit*
+                double result = 0;
+                for(int k = 0; k <= (long)n/4; k++){
+                    long factor = (n - 4*k);
+                    result += (double)mul_exp2_mod(4, factor, 8*k + 1)/(8*k+1);
+                    result -= (double)mul_exp2_mod(2, factor, 8*k + 4)/(8*k+4);
+                    result -= (double)mul_exp2_mod(1, factor, 8*k + 5)/(8*k+5);
+                    result -= (double)mul_exp2_mod(1, factor, 8*k + 6)/(8*k+6);
                 }
+                if (result < 0) result -= floor(result);
+                return (int)(result * 2);
+            }
+
+            // e^x = sum_{k=0}^\inf = x^0/0! + x^1/1! + x^2/2! + x^3/3! + x^3/6! + ...
 
-            };
+            //class exponential {
+            //    private:
+            //    // would be nice to interoperate between long, int, and boost::real::real,
+            //    // and have ctors from the integral types
+            //    boost::real::real k_prev = boost::real::real_explicit("0");
+            //    boost::real::real const * const x_ptr;
+            //    boost::real::real last_term; // x^kn / kn!
+            //    boost::real::real current_value; // summation from k0 to k_n, with precision digits
+
+            //    public:
+            //    exponential(boost::real::real &x) : x_ptr(&x) {
+            //        last_term = boost::real::real ("1");
+            //        current_value = boost::real::real ("1");
+            //    };
+
+            //    int get_nth_digit(unsigned int n) {
+            //        boost::real::real one = boost::real::real_explicit("1");
+            //        // if n < k_prev, reset
+
+            //        boost::real::real min_bound;
+            //        std::get<boost::real::real_explicit>(min_bound.get_real_number()).digits = {1};
+            //        std::get<boost::real::real_explicit>(min_bound.get_real_number()).exponent = 1-n;
+
+            //        // keep getting terms from the taylor series until the terms go below our precision bound
+            //        while((last_term > min_bound) || (last_term == min_bound)) {
+            //            last_term = last_term * (*x_ptr) / (k_prev + one);
+            //            current_value = current_value + last_term;
+            //        }
+
+            //        return 0;
+            //    }
+
+            //};
         }
     }
 }
diff --git a/test/irrational_helpers_test.cpp b/test/irrational_helpers_test.cpp
new file mode 100644
index 0000000..8375751
--- /dev/null
+++ b/test/irrational_helpers_test.cpp
@@ -0,0 +1,14 @@
+#include <catch2/catch.hpp>
+#include <real/real.hpp>
+#include <real/irrational_helpers.hpp>
+#include <string>
+
+TEST_CASE( "Compute irrational numbers", "[irrational]" ) {
+
+    SECTION("2000 bits of pi computed with BBP algorithm") {
+        std::string pi("1100100100001111110110101010001000100001011010001100001000110100110001001100011001100010100010111000000011011100000111001101000100101001000000100100111000001000100010100110011111001100011101000000001000001011101111101010011000111011000100111001101100100010010100010100101000001000011110011000111000110100000001001101110111101111100101010001100110110011110011010011101001000011000110110011000000101011000010100110110111110010010111110001010000110111010011111110000100110101011011010110110101010001110000100100010111100100100001011011010101110110011000100101111001111110110001101111010001001100010000101110100110100110001101111110110101101011000010111111111101011100101101101111010000000110101101111110110111101110001110000110101111111011010110101000100110011111101001011010111010011111001001000001000101111100010010110001111111100110010010010010100001100110010100011110110011100100010110110011110111000010000000000111110010111000101000010110001110111111000001011001100011011010010010000011011000011100010101011101001110011010011010010001011000111111101010001111110100100100110011110101111110000011011001010101110100100011110111001010001110101101100101100001110001100010111100110101011000100000100001010101001010111011100111101101010100101001000001110111000010010110100101100110110101100111000011000011010101001110010010101011110010011000000001001111000101110100011011000000100011001010000110000010000101111100001100101001000001011110010001100010111000110110110011100011101111100011100111100111011100101100000110000000111010000110000000111001101100100111100000111010001011101100000001111010001010001111101101011100010101011101111100000110111101001100010100101100100111011110001010111100101111110110100101010101100000010111000110000011100110010101010010010111110011101010100101010110101011100101000101011101001000100110000110001001100011111010000001010001000000010101011100101000111001011010100010101010101011000100001011011010110100110011000101110000110100000100010100000111101000110011101010000101010100");
+        for(int i = 1; i <= pi.length(); i++){
+            CHECK(boost::real::irrational::pi_binary_get_nth_digit(i) == (pi[i-1] - '0'));
+        }
+    }
+}