Skip to content

Commit c3b754f

Browse files
committed
Switched to MPIR for supporting larger decimal expansion on Windows.
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.
1 parent aac1538 commit c3b754f

File tree

19 files changed

+3163
-3795
lines changed

19 files changed

+3163
-3795
lines changed

ChudnovskyPiBS.cpp

Lines changed: 77 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -3,9 +3,13 @@
33
#include <barrier>
44
#include "ChudnovskyPiBS.h"
55
#include <sstream>
6+
#include "MPIR/gmp.h"
7+
#include "MPIR/gmpxx.h"
68

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

11+
#define HUNDRED_MILLION 100000000
12+
913
//Global variables for configuration:
1014
const int NUM_THREADS_IN_CPU = std::thread::hardware_concurrency();
1115
const int NUM_MAIN_WORKER_THREADS = NUM_THREADS_IN_CPU;
@@ -22,19 +26,51 @@ ChudnovskyPiBS::ChudnovskyPiBS(unsigned long _digits)
2226
throw _EXCEPTION_; //precision can't be handled by long double.
2327
N = (unsigned long)(digits / DIGITS_PER_TERM + 1);
2428

25-
mpz_ui_pow_ui(one_squared.get_mpz_t(), 10, 2 * digits);
26-
27-
//one_squared.get_str();
28-
29-
mpz_class thousandandfive__times__one_squared = 10005 * one_squared;
30-
mpz_sqrt(sqrtC.get_mpz_t(), thousandandfive__times__one_squared.get_mpz_t());
31-
32-
//sqrtC.get_str();
33-
3429
mpz_ui_pow_ui(intBigC3_OVER_24.get_mpz_t(), C, 3);
3530
mpz_class twentyfour = 24;
3631
mpz_fdiv_q(intBigC3_OVER_24.get_mpz_t(), intBigC3_OVER_24.get_mpz_t(), twentyfour.get_mpz_t());
3732
//intBigC3_OVER_24.get_str();
33+
if (digits <= HUNDRED_MILLION)
34+
{
35+
mpz_ui_pow_ui(one_squared.get_mpz_t(), 10, 2*_digits);
36+
mpz_class thousandandfive__times__one_squared = 10005 * one_squared;
37+
mpz_sqrt(sqrtC.get_mpz_t(), thousandandfive__times__one_squared.get_mpz_t());
38+
//std::cout << "Needed strlen:" << sqrtC.get_str().length() << std::endl;
39+
}
40+
else
41+
{
42+
mp_bitcnt_t prec = digits * log2l(10); //or is it 8 * (digits + 4)??
43+
SHIFTER.set_prec(prec);
44+
sqrtCF.set_prec(prec);
45+
46+
mpf_class ten = mpf_class(10);
47+
mpf_pow_ui(SHIFTER.get_mpf_t(), ten.get_mpf_t(), digits);
48+
//std::cout << "Prec:" << SHIFTER.get_prec() << std::endl;
49+
mpf_sqrt_ui(sqrtCF.get_mpf_t(),10005);
50+
//sqrtCF += 0.1 / SHIFTER;
51+
sqrtCF *= SHIFTER;
52+
mpf_trunc(sqrtCF.get_mpf_t(), sqrtCF.get_mpf_t());
53+
mp_exp_t exp = 0;
54+
std::string tempMiddleWare = sqrtCF.get_str(exp, 10, _digits + 4);
55+
//std::cout << "Gotten strlen:" << tempMiddleWare.length() << std::endl;
56+
//std::cout << "Decimal pt location:" << exp << std::endl;
57+
sqrtC = mpz_class(tempMiddleWare);
58+
//sqrtC.get_str();
59+
/*
60+
////ten = SHIFTER;
61+
////for (int i = 1; i < 10; i++) //this is basically doing (10^digits/10)^10 (by power rule, it becomes 10^digits)
62+
//// SHIFTER *= ten;
63+
////std::cout << one_squared.get_str() << std::endl;
64+
////one_squared *= one_squared; //this makes it 10^(digits+digits) = 10^(2*digits) which is what we wanted.
65+
////mpz_class lol;
66+
////mpz_ui_pow_ui(lol.get_mpz_t(), 10, 2*_digits);
67+
////if (lol != one_squared)
68+
////{
69+
//// std::cout << one_squared.get_str() << std::endl;
70+
//// std::cout << lol.get_str() << std::endl;
71+
////}
72+
*/
73+
}
3874
}
3975

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

272+
void ChudnovskyPiBS::getSqrtC(unsigned long digits)
273+
{
274+
if (digits <= HUNDRED_MILLION)
275+
{
276+
mpz_ui_pow_ui(one_squared.get_mpz_t(), 10, 2 * digits);
277+
}
278+
else
279+
{
280+
mpz_class ten = 10;
281+
mpz_pow_ui(one_squared.get_mpz_t(), ten.get_mpz_t(), digits / 10); //10^(digits/10)
282+
/*mpz_pow_ui(one_squared.get_mpz_t(), one_squared.get_mpz_t(), 10);*/
283+
ten = one_squared;
284+
for (int i = 1; i < 10; i++) //this is basically doing (10^digits/10)^10 (by power rule, it becomes 10^digits)
285+
one_squared *= ten;
286+
//std::cout << one_squared.get_str() << std::endl;
287+
one_squared *= one_squared; //this makes it 10^(digits+digits) = 10^(2*digits) which is what we wanted.
288+
//mpz_class lol;
289+
//mpz_ui_pow_ui(lol.get_mpz_t(), 10, 2*_digits);
290+
//if (lol != one_squared)
291+
//{
292+
// std::cout << one_squared.get_str() << std::endl;
293+
// std::cout << lol.get_str() << std::endl;
294+
//}
295+
}
296+
mpz_class thousandandfive__times__one_squared = 10005 * one_squared;
297+
mpz_sqrt(sqrtC.get_mpz_t(), thousandandfive__times__one_squared.get_mpz_t());
298+
}
299+
236300
int ChudnovskyPiBS::getTotalNumThreadsFromUsefulNumThreads(int usefulThreadCountWanted)
237301
{
238302
if (usefulThreadCountWanted == 2)
@@ -257,7 +321,8 @@ mpz_class ChudnovskyPiBS::calculatePi()
257321
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.
258322
//BSResult.Q.get_str();
259323
//BSResult.T.get_str();
260-
//return mpz_fdiv_q((Q * 426880 * sqrtC) / T
324+
325+
/*getSqrtC(digits);*/
261326
mpz_class result = (BSResult.Q * 426880 * sqrtC);
262327
mpz_fdiv_q(result.get_mpz_t(), result.get_mpz_t(), BSResult.T.get_mpz_t());
263328

@@ -266,6 +331,8 @@ mpz_class ChudnovskyPiBS::calculatePi()
266331
}
267332

268333

334+
335+
269336
/* DEPRECATED/OLD MULTITHREADING SOLUTION
270337
mpz_class ChudnovskyPiBS::calculatePi()
271338
{

ChudnovskyPiBS.h

Lines changed: 11 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,8 @@
33
#include <math.h>
44
#include <stdexcept>
55
#include <barrier>
6-
#include "GMP/gmpxx.h"
6+
#include "MPIR/gmp.h"
7+
#include "MPIR/gmpxx.h"
78

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

2526
unsigned long N;
2627
unsigned long digits;
27-
mpz_class one_squared;
28-
mpz_class sqrtC;
28+
mpz_class one_squared = mpz_class(-1);
29+
mpf_class SHIFTER = mpf_class(-1);
2930

30-
mpz_class oneClass = 1;
31-
mpz_ptr one = oneClass.get_mpz_t();
31+
mpz_class sqrtC = mpz_class(-1);
32+
mpf_class sqrtCF = mpf_class(-1);
33+
34+
mpz_class digitOneClass = mpz_class(1);
35+
mpz_ptr one = digitOneClass.get_mpz_t();
3236

3337
/// <summary>
3438
/// Computes the terms for binary splitting the Chudnovsky infinite series.
@@ -48,7 +52,7 @@ class ChudnovskyPiBS //Calculates Pi using Chudnovsky algorithm with Binary Spli
4852
void directlyCompute__P_Q_T__from_A_to_AplusOne(mpz_class& a, mpz_class& Pab, mpz_class& Qab, mpz_class& Tab);
4953
bsReturn bs_multithreaded(mpz_class a, mpz_class b, int threadCount);
5054
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.
51-
55+
void getSqrtC(unsigned long digits);
5256
public:
5357
/// <summary>
5458
/// Constructor for ChudnovskyPiBS Class.

GMP/DLLs/libgmp-10.dll

-642 KB
Binary file not shown.

GMP/DLLs/libgmpxx-4.dll

-30 KB
Binary file not shown.

GMP/Linker Files/libgmp.dll.a

-554 KB
Binary file not shown.

GMP/Linker Files/libgmpxx.dll.a

-70.2 KB
Binary file not shown.
-750 KB
Binary file not shown.
-62.5 KB
Binary file not shown.
-536 KB
Binary file not shown.
-71.5 KB
Binary file not shown.

0 commit comments

Comments
 (0)