Skip to content

Commit

Permalink
Switched to MPIR for supporting larger decimal expansion on Windows.
Browse files Browse the repository at this point in the history
Also, added supporting code in ChudnovskyPiBS class for what I thought would be the most efficient way to handle 10^9 decimal places at the time.
  • Loading branch information
Nomi committed Dec 3, 2022
1 parent aac1538 commit c3b754f
Show file tree
Hide file tree
Showing 19 changed files with 3,163 additions and 3,795 deletions.
87 changes: 77 additions & 10 deletions ChudnovskyPiBS.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -3,9 +3,13 @@
#include <barrier>
#include "ChudnovskyPiBS.h"
#include <sstream>
#include "MPIR/gmp.h"
#include "MPIR/gmpxx.h"

#include <Windows.h> //this should be the last import somehow?

#define HUNDRED_MILLION 100000000

//Global variables for configuration:
const int NUM_THREADS_IN_CPU = std::thread::hardware_concurrency();
const int NUM_MAIN_WORKER_THREADS = NUM_THREADS_IN_CPU;
Expand All @@ -22,19 +26,51 @@ ChudnovskyPiBS::ChudnovskyPiBS(unsigned long _digits)
throw _EXCEPTION_; //precision can't be handled by long double.
N = (unsigned long)(digits / DIGITS_PER_TERM + 1);

mpz_ui_pow_ui(one_squared.get_mpz_t(), 10, 2 * digits);

//one_squared.get_str();

mpz_class thousandandfive__times__one_squared = 10005 * one_squared;
mpz_sqrt(sqrtC.get_mpz_t(), thousandandfive__times__one_squared.get_mpz_t());

//sqrtC.get_str();

mpz_ui_pow_ui(intBigC3_OVER_24.get_mpz_t(), C, 3);
mpz_class twentyfour = 24;
mpz_fdiv_q(intBigC3_OVER_24.get_mpz_t(), intBigC3_OVER_24.get_mpz_t(), twentyfour.get_mpz_t());
//intBigC3_OVER_24.get_str();
if (digits <= HUNDRED_MILLION)
{
mpz_ui_pow_ui(one_squared.get_mpz_t(), 10, 2*_digits);
mpz_class thousandandfive__times__one_squared = 10005 * one_squared;
mpz_sqrt(sqrtC.get_mpz_t(), thousandandfive__times__one_squared.get_mpz_t());
//std::cout << "Needed strlen:" << sqrtC.get_str().length() << std::endl;
}
else
{
mp_bitcnt_t prec = digits * log2l(10); //or is it 8 * (digits + 4)??
SHIFTER.set_prec(prec);
sqrtCF.set_prec(prec);

mpf_class ten = mpf_class(10);
mpf_pow_ui(SHIFTER.get_mpf_t(), ten.get_mpf_t(), digits);
//std::cout << "Prec:" << SHIFTER.get_prec() << std::endl;
mpf_sqrt_ui(sqrtCF.get_mpf_t(),10005);
//sqrtCF += 0.1 / SHIFTER;
sqrtCF *= SHIFTER;
mpf_trunc(sqrtCF.get_mpf_t(), sqrtCF.get_mpf_t());
mp_exp_t exp = 0;
std::string tempMiddleWare = sqrtCF.get_str(exp, 10, _digits + 4);
//std::cout << "Gotten strlen:" << tempMiddleWare.length() << std::endl;
//std::cout << "Decimal pt location:" << exp << std::endl;
sqrtC = mpz_class(tempMiddleWare);
//sqrtC.get_str();
/*
////ten = SHIFTER;
////for (int i = 1; i < 10; i++) //this is basically doing (10^digits/10)^10 (by power rule, it becomes 10^digits)
//// SHIFTER *= ten;
////std::cout << one_squared.get_str() << std::endl;
////one_squared *= one_squared; //this makes it 10^(digits+digits) = 10^(2*digits) which is what we wanted.
////mpz_class lol;
////mpz_ui_pow_ui(lol.get_mpz_t(), 10, 2*_digits);
////if (lol != one_squared)
////{
//// std::cout << one_squared.get_str() << std::endl;
//// std::cout << lol.get_str() << std::endl;
////}
*/
}
}

bsReturn ChudnovskyPiBS::bs(mpz_class a, mpz_class b) //clearly thread safe because nothing from outside the function is written to.
Expand Down Expand Up @@ -233,6 +269,34 @@ bsReturn ChudnovskyPiBS::bs_multithreaded_barrier(mpz_class a, mpz_class b, int
return result;
}

void ChudnovskyPiBS::getSqrtC(unsigned long digits)
{
if (digits <= HUNDRED_MILLION)
{
mpz_ui_pow_ui(one_squared.get_mpz_t(), 10, 2 * digits);
}
else
{
mpz_class ten = 10;
mpz_pow_ui(one_squared.get_mpz_t(), ten.get_mpz_t(), digits / 10); //10^(digits/10)
/*mpz_pow_ui(one_squared.get_mpz_t(), one_squared.get_mpz_t(), 10);*/
ten = one_squared;
for (int i = 1; i < 10; i++) //this is basically doing (10^digits/10)^10 (by power rule, it becomes 10^digits)
one_squared *= ten;
//std::cout << one_squared.get_str() << std::endl;
one_squared *= one_squared; //this makes it 10^(digits+digits) = 10^(2*digits) which is what we wanted.
//mpz_class lol;
//mpz_ui_pow_ui(lol.get_mpz_t(), 10, 2*_digits);
//if (lol != one_squared)
//{
// std::cout << one_squared.get_str() << std::endl;
// std::cout << lol.get_str() << std::endl;
//}
}
mpz_class thousandandfive__times__one_squared = 10005 * one_squared;
mpz_sqrt(sqrtC.get_mpz_t(), thousandandfive__times__one_squared.get_mpz_t());
}

int ChudnovskyPiBS::getTotalNumThreadsFromUsefulNumThreads(int usefulThreadCountWanted)
{
if (usefulThreadCountWanted == 2)
Expand All @@ -257,7 +321,8 @@ mpz_class ChudnovskyPiBS::calculatePi()
bsReturn BSResult = bs_multithreaded(0, N, TOTAL_THREAD_COUNT); //bs_multithreaded(0, N, TOTAL_THREAD_COUNT); //bs_multithreaded_barrier(0, N, TOTAL_THREAD_COUNT, 0); //bs(0, N); //apparently Q and T gotten are wrong.
//BSResult.Q.get_str();
//BSResult.T.get_str();
//return mpz_fdiv_q((Q * 426880 * sqrtC) / T

/*getSqrtC(digits);*/
mpz_class result = (BSResult.Q * 426880 * sqrtC);
mpz_fdiv_q(result.get_mpz_t(), result.get_mpz_t(), BSResult.T.get_mpz_t());

Expand All @@ -266,6 +331,8 @@ mpz_class ChudnovskyPiBS::calculatePi()
}




/* DEPRECATED/OLD MULTITHREADING SOLUTION
mpz_class ChudnovskyPiBS::calculatePi()
{
Expand Down
18 changes: 11 additions & 7 deletions ChudnovskyPiBS.h
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,8 @@
#include <math.h>
#include <stdexcept>
#include <barrier>
#include "GMP/gmpxx.h"
#include "MPIR/gmp.h"
#include "MPIR/gmpxx.h"

//Helper structs:
struct bsReturn
Expand All @@ -19,16 +20,19 @@ class ChudnovskyPiBS //Calculates Pi using Chudnovsky algorithm with Binary Spli
private:
const long double C = 640320;
const long double C3_OVER_24 = powl(C,3)/24;
mpz_class intBigC3_OVER_24 = 0;
mpz_class intBigC3_OVER_24 = mpz_class(0);
const long double DIGITS_PER_TERM = log10l(C3_OVER_24 / 6 / 2 / 6);

unsigned long N;
unsigned long digits;
mpz_class one_squared;
mpz_class sqrtC;
mpz_class one_squared = mpz_class(-1);
mpf_class SHIFTER = mpf_class(-1);

mpz_class oneClass = 1;
mpz_ptr one = oneClass.get_mpz_t();
mpz_class sqrtC = mpz_class(-1);
mpf_class sqrtCF = mpf_class(-1);

mpz_class digitOneClass = mpz_class(1);
mpz_ptr one = digitOneClass.get_mpz_t();

/// <summary>
/// Computes the terms for binary splitting the Chudnovsky infinite series.
Expand All @@ -48,7 +52,7 @@ class ChudnovskyPiBS //Calculates Pi using Chudnovsky algorithm with Binary Spli
void directlyCompute__P_Q_T__from_A_to_AplusOne(mpz_class& a, mpz_class& Pab, mpz_class& Qab, mpz_class& Tab);
bsReturn bs_multithreaded(mpz_class a, mpz_class b, int threadCount);
bsReturn bs_multithreaded_barrier(mpz_class a, mpz_class b, int threadCount, int depth); //uses a barrier to wait for all main worker threads to spawn.

void getSqrtC(unsigned long digits);
public:
/// <summary>
/// Constructor for ChudnovskyPiBS Class.
Expand Down
Binary file removed GMP/DLLs/libgmp-10.dll
Binary file not shown.
Binary file removed GMP/DLLs/libgmpxx-4.dll
Binary file not shown.
Binary file removed GMP/Linker Files/libgmp.dll.a
Binary file not shown.
Binary file removed GMP/Linker Files/libgmpxx.dll.a
Binary file not shown.
Binary file removed GMP/Old_x86/CopyToOutputBinFolder/libgmp-10.dll
Binary file not shown.
Binary file removed GMP/Old_x86/CopyToOutputBinFolder/libgmpxx-4.dll
Binary file not shown.
Binary file removed GMP/Old_x86/Linker Files/libgmp.dll.a
Binary file not shown.
Binary file removed GMP/Old_x86/Linker Files/libgmpxx.dll.a
Binary file not shown.
Binary file added MPIR/DLLs/mpir.dll
Binary file not shown.
Binary file added MPIR/Linker files/mpir.lib
Binary file not shown.
Loading

0 comments on commit c3b754f

Please sign in to comment.