Skip to content

Implementation of pi_binary_get_nth_digit using BBP algorithm #47

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
firmart wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Implementation of pi_binary_get_nth_digit using BBP algorithm #47

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
e3e37e6
Add pi_binary_get_nth_digit() using BBP algorithm.
firmart Mar 30, 2020
2088794
add irrational_helpers test (including pi).
firmart Mar 30, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
99 changes: 61 additions & 38 deletions include/real/irrational_helpers.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
// }

//};
}
}
}
Expand Down
14 changes: 14 additions & 0 deletions test/irrational_helpers_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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'));
}
}
}