eTo compute mathematical constants and functions with high precision, we must express them using the fundamental arithmetic operations: addition, subtraction, multiplication, and division. This is usually achieved by finding a Taylor series expansion of the constant or function centered at a specific, well-understood point.
Taylor Series
A Taylor series is an infinite sum of terms that, when fully expanded, approximates a given function. When this sum approaches the actual function value, we say the series converges to the function. The Taylor series is developed around a specific point within the function’s domain and is used to estimate function values near that point. If the expansion is centered at 0, it is called a Maclaurin series. To construct a Taylor series, we must know the derivatives of the function. Given this information, the Taylor series centered at a point \(a\) is defined as:
\[ f(x) = f(a) + f'(a)(x - a) + \frac{f''(a)}{2!}(x - a)^2 + \cdots = \sum_{n=0}^{\infty} \frac{f^{(n)}(a)}{n!} (x - a)^n\]
The factorial of \(n\) is written as \(n!\). The expressions \(f(a)\), \(f'(a)\) and \(f''(a)\) represent the value of the function and its first and second derivatives, all evaluated at a fixed point \(a\). More generally, \(f^{(n)}(a)\) denotes the nth derivative of \(f\) evaluated at \(a\). By convention, the zero-order derivative of \(f\) is defined as \(f\) itself, and \((x − a)^0\) and \(0!\) are both defined to be 1.
Let’s walk through the example of constructing the Taylor series for \(e^x\). This is a straightforward case since all derivatives of \(e^x\) are all \(e^x\).
\[ e^x = \sum_{n=0}^{\infty} \frac{e^x}{n!} (x - a)^n\]
If we take the fixed point \(a\) as being 0 we get the Maclauren series for \(e^x\).
\[ e^x = \sum_{n=0}^{\infty} \frac{x^n}{n!} = 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} + \frac{x^4}{4!} + \cdots \]
We can create a simple C++ program to calculate the value of \(e^x\) using its Maclauren series expansion.
maclauren_exp.cpp
/**
* @file maclauren_exp.cpp
* @brief Program to compute \f(e^x\f).
*
* This program can compute \f(e^x\f) using the Maclauren series.
* \f[e^x = \sum_{n=0}^{\infty} \frac{x^n}{n!} = 1 + x + \frac{x^2}{2!}
* + \frac{x^3}{3!} + \frac{x^4}{4!} + \cdots\f]
*
* Syntax: maclauren_exp [x]
*
* @author Richard Lesh
* @date 2025-07-27
* @version 1.0
*/
#include <chrono>
#include <cmath>
#include <iostream>
using namespace std;
/**
* @brief Computes the factorial of a non-negative integer.
*
* This function recursively computes the factorial \f(n!\f) of a given
* non-negative integer `n`. The factorial of 0 and 1 is defined as 1. For
* values greater than 1, the function returns `n * factorial(n - 1)`.
*
* @param n A non-negative integer whose factorial is to be computed. Max. 20.
* @return The factorial of the input number `n`. Returns 1 if `n` is 0 or 1.
*
* @note The function does not perform overflow checking. Results are undefined
* for values of `n` that cause overflow in the return type `long`, i.e.,
* n > 20.
*/
// NOLINTNEXTLINE(misc-no-recursion)
long factorial(long n) {
if (n <= 1) return 1;
return n * factorial(n - 1);
}
/**
* @brief Approximates the exponential function \f(e^x\f) using a Maclaurin series.
*
* This function computes an approximation of \f(e^x\f) by evaluating the Maclaurin
* series expansion up to 20 terms or until the contribution of a term becomes
* insignificant relative to the accumulated sum (below a relative tolerance of
* \f(10^{-16}\f)).
*
* The stopping condition ensures the result does not include terms that are
* insignificant due to the limits of double-precision floating-point accuracy
* (approximately 16 decimal digits). But due to the limitations of double
* precision floating point range when computing the factorial, at maximum 20
* terms can be computed resulting in errors in the approximation for values of
* x not near 0.
*
* @param x The exponent value for which \f(e^x\f) is to be approximated.
* @return An approximation of \f(e^x\f) computed using the Maclaurin series.
*
* @note The function prints the number of terms used before early termination.
*/
double maclauren_exp(double x) {
double sum{0};
// Limit terms to i <= 20 because of max. factorial
for (int i = 0; i <= 20; ++i) {
double term = pow(x, i) / factorial(i);
sum += term;
// Break if the relative error (estimated to be the ratio of the
// term / sum) is less than tolerance, \f(10^{-16}\f). We use this tolerance
// because doubles have about 16 digits of accuracy.
if (abs(term) <= sum * 1e-16) {
cout << "Terms: " << i << endl;
break;
}
}
return sum;
}
/**
* @brief Main function.
*
* The entry point of the program.
*
* Syntax: maclauren_exp [x]
*
* @return int Returns 0 on successful execution.
*/
int main(int argc, char **argv) {
double x{0.5};
if (argc > 1) {
x = atof(argv[1]);
}
// Use the high_resolution_clock to time execution
auto start = chrono::high_resolution_clock::now();
double exp = maclauren_exp(x);
auto stop = chrono::high_resolution_clock::now();
chrono::duration<double> elapsed = stop - start;
printf("x = %.16g\n", x);
printf("e^x = %.16g\n", exp);
printf("std::exp(x) = %.16g\n", std::exp(x));
// We can compute the error in our approximation using the standard library
// function exp(x).
printf("error = %g%%\n", abs(exp - std::exp(x)) / std::exp(x) * 100);
printf("Computation Time: %g seconds\n", elapsed.count());
return 0;
}
% maclauren_exp -1
Terms: 19
x = -1
e^x = 0.3678794411714424
std::exp(x) = 0.3678794411714423
error = 3.0179e-14%
Computation Time: 1.275e-05 seconds
% maclauren_exp 0.5
x = 0.5
Terms: 15
e^x = 1.648721270700128
std::exp(x) = 1.648721270700128
error = 2.69354e-14%
Computation Time: 4.459e-06 seconds
% maclauren_exp 1
x = 1
Terms: 18
e^x = 2.718281828459046
std::exp(x) = 2.718281828459045
error = 1.63371e-14%
Computation Time: 4.917e-06 seconds
% maclauren_exp 1.4
x = 1.4
Terms: 20
e^x = 4.055199966844675
std::exp(x) = 4.055199966844675
error = 0%
Computation Time: 6.542e-06 seconds
% maclauren_exp 5
x = 5
e^x = 148.4131470673818
std::exp(x) = 148.4131591025766
error = 8.10925e-06%
Computation Time: 1.084e-06 seconds
% maclauren_exp 10
x = 10
e^x = 21991.48202566506
std::exp(x) = 22026.46579480672
error = 0.158826%
Computation Time: 1.25e-06 seconds
This program’s accuracy is limited to about 16 decimal digits due to its reliance on IEEE 754 double-precision floating-point numbers. Additionally, it can use at most 21 terms, since \(20!\) is the largest factorial that can be accurately represented within the bounds of IEEE 754 doubles.
When we run the computation with values near 0 — the point around which the series is expanded — we obtain accurate results using only a few terms. However, as we move farther from \(x = 0\), more terms are needed to maintain the same level of accuracy. By the time we reach \(x = 1.4\), we’ve hit the limit of what can be computed accurately with the 20-term cap we’ve imposed, due to the constraints of the factorial function. For values like 5 and 10, the percentage error continues to grow.
With the basic version of the \(exp(x)\) function working using double precision, we can now extend the program to support arbitrary-precision arithmetic using the GNU Multiple Precision Arithmetic Library (GMP). This library allows computations with as many digits of precision as needed.
By using GMP’s mpz_class for arbitrary-size integers, we can implement a factorial function without the 20-term input limitation. And by using the mpreal class for high-precision floating-point numbers, we can create a maclauren_exp() function that computes \(e^x\) to any desired level of accuracy.
maclauren_exp_mp.cpp
/**
* @file maclauren_exp.cpp
* @brief Program to compute \f(e^x\f).
*
* This program can compute \f(e^x\f) using the Maclauren series.
* \f[e^x = \sum_{n=0}^{\infty} \frac{x^n}{n!} = 1 + x + \frac{x^2}{2!}
* + \frac{x^3}{3!} + \frac{x^4}{4!} + \cdots\f]
*
* Syntax: maclauren_exp [x]
*
* @author Richard Lesh
* @date 2025-07-27
* @version 1.0
*/
#include <chrono>
#include <cmath>
#include <iostream>
using namespace std;
/**
* @brief Computes the factorial of a non-negative integer.
*
* This function recursively computes the factorial \f(n!\f) of a given
* non-negative integer `n`. The factorial of 0 and 1 is defined as 1. For
* values greater than 1, the function returns `n * factorial(n - 1)`.
*
* @param n A non-negative integer whose factorial is to be computed. Max. 20.
* @return The factorial of the input number `n`. Returns 1 if `n` is 0 or 1.
*
* @note The function does not perform overflow checking. Results are undefined
* for values of `n` that cause overflow in the return type `long`, i.e.,
* n > 20.
*/
// NOLINTNEXTLINE(misc-no-recursion)
long factorial(long n) {
if (n <= 1) return 1;
return n * factorial(n - 1);
}
/**
* @brief Approximates the exponential function \f(e^x\f) using a Maclaurin series.
*
* This function computes an approximation of \f(e^x\f) by evaluating the Maclaurin
* series expansion up to 20 terms or until the contribution of a term becomes
* insignificant relative to the accumulated sum (below a relative tolerance of
* \f(10^{-16}\f)).
*
* The stopping condition ensures the result does not include terms that are
* insignificant due to the limits of double-precision floating-point accuracy
* (approximately 16 decimal digits). But due to the limitations of double
* precision floating point range when computing the factorial, at maximum 20
* terms can be computed resulting in errors in the approximation for values of
* x not near 0.
*
* @param x The exponent value for which \f(e^x\f) is to be approximated.
* @return An approximation of \f(e^x\f) computed using the Maclaurin series.
*
* @note The function prints the number of terms used before early termination.
*/
double maclauren_exp(double x) {
double sum{0};
// Limit terms to i <= 20 because of max. factorial
for (int i = 0; i <= 20; ++i) {
double term = pow(x, i) / factorial(i);
sum += term;
// Break if the relative error (estimated to be the ratio of the
// term / sum) is less than tolerance, \f(10^{-16}\f). We use this tolerance
// because doubles have about 16 digits of accuracy.
if (abs(term) <= sum * 1e-16) {
cout << "Terms: " << i << endl;
break;
}
}
return sum;
}
/**
* @brief Main function.
*
* The entry point of the program.
*
* Syntax: maclauren_exp [x]
*
* @return int Returns 0 on successful execution.
*/
int main(int argc, char **argv) {
double x{0.5};
if (argc > 1) {
x = atof(argv[1]);
}
// Use the high_resolution_clock to time execution
auto start = chrono::high_resolution_clock::now();
double exp = maclauren_exp(x);
auto stop = chrono::high_resolution_clock::now();
chrono::duration<double> elapsed = stop - start;
printf("x = %.16g\n", x);
printf("e^x = %.16g\n", exp);
printf("std::exp(x) = %.16g\n", std::exp(x));
// We can compute the error in our approximation using the standard library
// function exp(x).
printf("error = %g%%\n", abs(exp - std::exp(x)) / std::exp(x) * 100);
printf("Computation Time: %g seconds\n", elapsed.count());
return 0;
}
% maclauren_exp_mp 1 16
x = 1
Terms: 18
e^x = 2.718281828459045
std::exp(x) = 2.718281828459045
Computation Time: 2.3666e-05 seconds
% maclauren_exp_mp 10 16
x = 10
Terms: 46
e^x = 22026.46579480672
std::exp(x) = 22026.46579480672
Computation Time: 3.5334e-05 seconds
% maclauren_exp_mp 1 50
x = 1
Terms: 42
e^x = 2.7182818284590452353602874713526624977572470937
std::exp(x) = 2.718281828459045
Computation Time: 4.2708e-05 seconds
% maclauren_exp_mp 1 1000
x = 1
Terms: 450
e^x = 2.718281828459045235360287471352662497757247093699959574966967627724076630
35354759457138217852516642742746639193200305992181741359662904357290033429526059
56307381323286279434907632338298807531952510190115738341879307021540891499348841
67509244761460668082264800168477411853742345442437107539077744992069551702761838
60626133138458300075204493382656029760673711320070932870912744374704723069697720
93101416928368190255151086574637721112523897844250569536967707854499699679468644
54905987931636889230098793127736178215424999229576351482208269895193668033182528
86939849646510582093923982948879332036250944311730123819706841614039701983767932
06832823764648042953118023287825098194558153017567173613320698112509961818815930
41690351598888519345807273866738589422879228499892086805825749279610484198444363
46324496848756023362482704197862320900216099023530436994184914631409343173814364
05462531520961836908887070167683964243781405927145635490613031072085103837505101
15747704171898610687396965521267154688957035035
std::exp(x) = 2.718281828459045
Computation Time: 0.00309975 seconds
At 16 decimal digits of precision, we can compute \(e^1\) and obtain the same result we previously achieved using double precision with 18 terms. When computing \(e^{10}\), we now get an accurate result because we’re no longer restricted to just 20 terms by the limitations of the factorial function. In fact, we can run the program to compute \(e^1\) with 50 or more digits of accuracy. However, as the number of desired digits increases, so does the number of terms required, which in turn increases computation time. For example, on a 2023 MacBook Pro M2, it takes approximately 33.24 seconds to compute 50,000 digits of \(e^1\).
There are several ways to optimize this process. The first improvement comes from recognizing that we don’t need to recompute the entire factorial from scratch for each term. Instead, we can reuse the factorial from the previous term and multiply it by the next integer. This optimization reduces the runtime for computing 50,000 digits of \(e^1)\) to 14.76 seconds.
Another potential optimization is to replace the standard power function with a custom integer exponentiation function using exponentiation by squaring, since the exponent in each term is an integer. After implementing this ipow() function and testing it, we find that the performance remains unchanged — indicating that mpreal::pow() already uses exponentiation by squaring internally.
This highlights a fundamental principle of optimization: avoid premature optimization. Focus first on writing correct and working code. Once that’s done, measure performance, apply one optimization at a time, and verify that each change actually improves efficiency.
We can apply a final optimization, similar to the one used for the factorial. Observing the numerator in each term, we see that it increases by just one factor of \(x\) compared to the previous term. Rather than recalculating the power from scratch with pow(), we can carry the numerator forward by multiplying it by \(x\) on each iteration. This reduces the numerator computation to a single multiplication per loop, in contrast to the multiple operations required by the ipow() function. Implementing this change yields a modest performance improvement, reducing the computation time to 13.38 seconds.
maclauren_exp_mp_opt.cpp
/**
* @file maclauren_exp_mp_opt.cpp
* @brief Program to compute \f(e^x\f).
*
* This program can compute \f(e^x\f) using the Maclauren series with the
* GNU Multiple Precision (GMP) library. This optimized version uses the
* relationship between terms \f(t_n = \frac{x}{n} t_{n-1}\f) to reduce
* significantly the number of multiplication ops required to compute each term.
*
* \f[e^x = \sum_{n=0}^{\infty} \frac{x^n}{n!} = 1 + x + \frac{x^2}{2!}
* + \frac{x^3}{3!} + \frac{x^4}{4!} + \cdots\f]
*
* Syntax: maclauren_exp_mp_opt [x] [num_digits]
*
* get mpreal.h at https://github.com/advanpix/mpreal.git
*
* Compile with command:\n
* macOS:\n
* install homebrew https://brew.sh\n
* put mpreal.h in /opt/homebrew/include\n
* % brew install gcc\n
* % brew install gmp\n
* % brew install mpfr\n
* % g++ -std=c++20 -O3 -I/opt/homebrew/include -L/opt/homebrew/lib -o maclauren_exp_mp_opt maclauren_exp_mp_opt.cpp -lmpfr -lgmp
*
* Linux:\n
* put mpreal.h in /urs/local/include\n
* % sudo apt install libgmp-dev\n
* % sudo apt install libmpfr-dev\n
* % g++ -std=c++20 -O3 -o maclauren_exp_mp_opt maclauren_exp_mp_opt.cpp -lmpfr -lgmp
* @author Richard Lesh
* @date 2025-07-27
* @version 1.0
*/
#include <chrono>
#include <cmath>
#include <gmpxx.h> // GMP C++ Wrapper
#include <iomanip>
#include <iostream>
#include <mpreal.h> // MPFR C++ Wrapper
using namespace mpfr;
using namespace std;
long DIGITS = 16; // in decimal digits
long PRECISION; // in bits
using namespace std;
/**
* @brief Approximates the exponential function \f(e^x\f) using a Maclaurin series.
* The computation is carried out using a specified number of decimal DIGITS.
*
* This function computes an approximation of \f(e^x\f) by evaluating the Maclaurin
* series expansion until the contribution of a term becomes insignificant relative
* to the accumulated sum (below a relative tolerance of \f(10^{-DIGITS}\f)).
*
* Uses the relationship between terms \f(t_n = \frac{x}{n} t_{n-1}\f) to reduce
* significantly the number of multiplication ops required to compute each term.
*
* @param x The exponent value for which \f(e^x\f) is to be approximated.
* @return An approximation of \f(e^x\f) computed using the Maclaurin series.
*
* @note The function prints the number of terms used to compute the approximation.
*/
mpreal maclauren_exp(const mpreal& x) {
mpreal sum{1};
mpreal tolerance = pow(10.0, -DIGITS);
mpz_class factorial_mpz{1};
mpreal power = 1;
for (int i = 1; i < numeric_limits<int>::max(); ++i) {
// new term is the previous term times \f(\frac{x}{i}\f)
factorial_mpz *= i;
power *= x;
mpreal term = power / factorial_mpz.get_mpz_t();
sum += term;
// Break if the relative error (estimated to be the ratio of the
// term / sum) is less than tolerance, \f(10^{-DIGITS}\f).
if (abs(term) <= sum * tolerance) {
cout << "Terms: " << i << endl;
break;
}
}
return sum;
}
/**
* @brief Main function.
*
* The entry point of the program.
*
* Syntax: maclauren_exp_mp_opt [x] [num_digits]
*
* @return int Returns 0 on successful execution.
*/
int main(int argc, char **argv) {
if (argc > 2) {
// Using double conversion so num_digits can be in exponential notation
// such as 1e3 or 1e6
DIGITS = static_cast<long>(atof(argv[2]));
}
// Precision in bits (~3.3 bits per decimal digit)
PRECISION = static_cast<long>(DIGITS * 3.32193 + 16);
mpreal::set_default_prec(PRECISION); // Set default precision in bits
mpreal x{1.0};
if (argc > 1) {
x = mpreal(argv[1]);
}
// Use the high_resolution_clock to time execution
auto start = chrono::high_resolution_clock::now();
mpreal exp = maclauren_exp(x);
auto stop = chrono::high_resolution_clock::now();
chrono::duration<double> elapsed = stop - start;
cout << "x = " << x << endl;
cout << "e^x = " << exp.toString(DIGITS) << endl;
cout << "std::exp(x) = " << setprecision(16) << std::exp(x.toDouble()) << endl;
cout << "Computation Time: " << setprecision(3) << elapsed.count() << " seconds" << endl;
return 0;
}
% maclauren_exp_mp_opt 1 50000
x = 1
Terms: 13525
e^x = 2.718281828459045235360287471352662497757247093699959574966967627724076630
35354759457138217852516642742746639193200305992181741359662904357290033429526059
56307381323286279434907632338298807531952510190115738341879307021540891499348841
...
14492534279922013463805383423696437674288625951401461782018107341005654667082368
54312816339049676558789901487477972479202502227218169405159042170892104287552188
65830860845270842392865259753614629003778016700165467168160534329290757303146656
248581
std::exp(x) = 2.718281828459045
Computation Time: 13.383637625 seconds
Binary Splitting
Another way to speed up the computation is to find an algorithm that avoids division as much as possible. Of the four arithmetic operations, division is by far the most costly to compute. One such optimization technique is known as Binary Splitting.
Binary splitting is based on recursively breaking down a sum like:
\[ S = \sum_{k=0}^{n} \frac{a_k}{b_k} \]
into partial sums over smaller ranges, then combining them using numerators and denominators separately. It’s most efficient when:
- Each term is a rational number
- The numerator and denominator can be computed recursively
- You can express the entire sum as a single rational number \(\frac{P}{Q}\)
We start with a series in summation notation that represents a partial sum over the half open interval \([i,j)\):
\[ S_{(i,j)} = \sum_{n=i}^{j-1} \frac{a_n}{b_n} = \frac{a_i}{b_i} + \frac{a_{i+1}}{b_{i+1}} + \cdots + \frac{a_{j-1}}{b_{j-1}} \]
For the simple case of \(j - i = 1\) we have...
\[ S_{(i,i+1)} = \frac{a_i}{b_i}\]
For the case when \(j - i = 2\) we can convert the two terms to a single term over a common denominator.
\[ S_{(i,i+2)} = \frac{a_i}{b_i} + \frac{a_{i+1}}{b_{i+1}} = \frac{a_ib_{i+1} + a_{i+1}b_i}{b_ib_{i+1}}\]
The pseudocode for our splitting algorithm then becomes:
split(a, b):
if b - a == 1:
return (a_i, b_i)
else:
m = (a + b) / 2
(P1, Q1) = Split(a, m)
(P2, Q2) = Split(m, b)
return (P1 * Q2 + P2 * Q1, Q1 * Q2)
So if we want to compute any partial sum we simply call split(a,b) and it recursively subdivides the problem into simpler problems. To compute the first \(n\) terms of our series we have:
\[ (P, Q) = split(0,n)\] \[S_{(0,n)} = \frac{P}{Q} \]
Compared to term-by-term summation:
- Avoids repeated divisions
- Reduces total number of operations from \(O(n^2)\) to \(O(n \log n)\)
- Uses integer arithmetic (exact) and delays floating-point conversion
One key difference when using binary splitting is that we must determine in advance how many terms to compute. Unlike the sequential approach to evaluating a Taylor series — where we can stop once the terms become smaller than a specified tolerance — binary splitting requires a fixed number of terms upfront. For our function, this means:
\[f(x) = P_n(x) + R_n(x)\]
where
\[P_n(x) = \sum_{k=0}^{n} \frac{f^{(k)}(a)}{k!}(x - a)^k\]
is the Taylor polynomial of degree n, and
\[R_n(x) = \frac{f^{(n+1)}(\xi)}{(n+1)!}(x - a)^{n+1}\]
Here \(R_n(x)\) is the Lagrange form of the remainder. The Lagrange form of the remainder expresses the error in a Taylor series approximation and involves some value \(\xi\) within the open interval between \(a\) and \(x\). While \(\xi\) is not known precisely, its existence is guaranteed by the Mean Value Theorem. In the case of the Maclaurin series for \(e^x\), we have \(a = 0\).
Since the \(n+1\)st derivative of \(e^x\) is \(e^x\), the maximum value of \(f^{(n+1)}(\xi)\) occurs at \(\xi = x\) for \(x>0\), and at \(\xi = 0\) for \(x<=0\). This is because \(e^x\) is strictly increasing, so its maximum on the interval \([0, x]\) (for \(x > 0\)) is \(e^x\), and on the interval \([x, 0]\) (for \(x < 0\)) is \(e^0 = 1\).
This observation allows us to establish an upper bound on the error of the degree-n Taylor polynomial for \(e^x\) as:
\[R_n(x) = \frac{\max(1, e^x)}{(n+1)!}(x)^{n+1}\]
If we want the Taylor polynomial to be accurate to \(D\) significant decimal digits we write the relative error as:
\[ \left| \frac{R_n(x)}{e^x} \right| = \left| \frac{\max(1, e^x)}{(n+1)!}\frac{(x)^{n+1}}{e^x} \right| <= 10^{-D}\]
\[ \log_{10} \left| \frac{\max(1, e^x)}{(n+1)!}\frac{(x)^{n+1}}{e^x} \right| <= -D \]
\[ \log_{10}(max(1,e^x)) + (n+1)\log_{10}(\left| x \right|) - log_{10}((n+1)!) - x\log_{10}(e) <= -D \]
Using Stirling's approximation for the log of large factorials (\(i > 10)\):
\[ \log_{10}(i!) \approx i \cdot \log_{10}(\frac{i}{e}) + \frac{1}{2} \log_{10}(2\pi i)\]
and simplifying we have:
\[ max(0,x \cdot log_{10}(e)) + (n+1)\log_{10}(\left| x \right|) - (n + 1) ( \log_{10}(\frac{n + 1}{e}) - \frac{1}{2} \log_{10}(2\pi (n + 1)) - x\log_{10}(e) <= -D \]
While this looks complicated, we can easily solve for \(n\) given \(D\) and \(x\) using a binary search algorithm to find the zero of the function:
\[ f(n, D, x) = -max(0,x \cdot log_{10}(e)) - (n+1)\log_{10}(\left| x \right|) + (n + 1) ( \log_{10}(\frac{n + 1}{e}) + \frac{1}{2} \log_{10}(2\pi (n + 1)) + x\log_{10}(e) - D \]
The pseudocode for the binary search is:
search_for_zero(D, x, low = 20, high = 1e9):
while high - low > 1:
mid = (low + high) / 2
y = f(mid, D, x)
if y > 0:
high = mid
else:
low = mid
end while
return high
Now we have all the parts to implement our \(e^x\) computation using binary splitting.
maclauren_exp_bs.cpp
/**
* @file maclauren_exp_bs.cpp
* @brief Program to compute \f(e^x\f).
*
* This program can compute \f(e^x\f) using the Maclauren series with the
* GNU Multiple Precision (GMP) library.
* It uses a binary splitting algorithm to potentially speed up the computation.
*
* Syntax: maclauren_exp_bs [x] [num_digits]
* x can be a fraction like 2/3 or in decimal notation 1.5
*
* get mpreal.h at https://github.com/advanpix/mpreal.git
*
* Compile with command:\n
* macOS:\n
* install homebrew https://brew.sh\n
* put mpreal.h in /opt/homebrew/include\n
* % brew install gcc\n
* % brew install gmp\n
* % brew install mpfr\n
* % g++ -std=c++20 -O3 -I/opt/homebrew/include -L/opt/homebrew/lib -o maclauren_exp_bs maclauren_exp_bs.cpp -lmpfr -lgmp
*
* Linux:\n
* put mpreal.h in /urs/local/include\n
* % sudo apt install libgmp-dev\n
* % sudo apt install libmpfr-dev\n
* % g++ -std=c++20 -O3 -o maclauren_exp_bs maclauren_exp_bs.cpp -lmpfr -lgmp
* @author Richard Lesh
* @date 2025-07-27
* @version 1.0
*/
#include <chrono>
#include <cmath>
#include <gmpxx.h> // GMP C++ Wrapper
#include <iomanip>
#include <iostream>
#include <mpreal.h> // MPFR C++ Wrapper
using namespace mpfr;
using namespace std;
long DIGITS = 16; // in decimal digits
long PRECISION; // in bits
using namespace std;
/**
* @brief Computes the factorial of a non-negative integer.
*
* This function recursively computes the factorial \f(n!\f) of a given
* non-negative integer `n`. The factorial of 0 and 1 is defined as 1. For
* values greater than 1, the function returns `n * factorial(n - 1)`.
*
* @param n A non-negative integer whose factorial is to be computed.
* @return The factorial of the input number `n`. Returns 1 if `n` is 0 or 1.
*
* @note The function does not overflow because arbitrary-precision integers are used.
*/
// NOLINTNEXTLINE(misc-no-recursion)
mpz_class factorial(long n) {
if (n <= 1) return 1;
return n * factorial(n - 1);
}
/**
* @brief Recursive binary splitting algorithm for computing a rational sum.
*
* This function recursively computes a pair of integers (P, Q) such that
* \f[
* \frac{P}{Q} = \sum_{k=a}^{b-1} \frac{x^k}{k!}
* \f]
* using a binary splitting strategy for improved performance and numerical stability.
*
* @param a Lower bound of the summation interval (inclusive).
* @param b Upper bound of the summation interval (exclusive).
* @param x Rational base of the power series (as an `mpq_class`).
* @param[out] P Resulting numerator of the partial sum, accumulated via binary splitting.
* @param[out] Q Resulting denominator of the partial sum, accumulated via binary splitting.
*
* @note The result \f$ \frac{P}{Q} \f$ represents the sum of the terms
* \f$ \frac{x^k}{k!} \f$ for \f$ k \in [a, b) \f$.
*
* @warning Assumes that `factorial(a)` is computable for the given range.
* May not perform well for very large ranges without optimizations.
*/
// NOLINTNEXTLINE(misc-no-recursion)
void binary_split(int a, int b, const mpq_class& x, mpz_class& P, mpz_class& Q) {
if (b - a == 1) {
mpz_pow_ui(P.get_mpz_t(), x.get_num().get_mpz_t(), a);
mpz_pow_ui(Q.get_mpz_t(), x.get_den().get_mpz_t(), a);
Q *= factorial(a);
} else {
int m = (a + b) / 2;
mpz_class P1, P2;
mpz_class Q1, Q2;
binary_split(a, m, x, P1, Q1);
binary_split(m, b, x, P2, Q2);
P = P1 * Q2 + P2 * Q1;
Q = Q1 * Q2;
}
}
/**
* @brief Approximates the number of digits contributed by a term in a power
* series expansion.
*
* This function models the relationship between the number of terms n,
* desired decimal digits of precision digits, and the evaluation point x
* using Stirling’s approximation. It returns the difference between the estimated
* number of digits the partial sum contributes and the required number of digits.
*
* The approximation is based on:
* \f[
* f(n, D, x) = -max(0,x \cdot log_{10}(e)) - (n+1)\log_{10}(\left| x \right|)
* + (n + 1) ( \log_{10}(\frac{n + 1}{e}) + \frac{1}{2} \log_{10}(2\pi (n + 1))
* + x\log_{10}(e) - D
* \f]
*
* @param n The number of terms used in the approximation.
* @param digits The target number of decimal digits for precision.
* @param x The evaluation point of the series.
* @return Estimated log10 error contribution; positive means more digits
* than needed.
*
* @note This function is intended to be used as a scoring or objective function
* in root-finding or binary search for optimal n.
*/
double digits_to_terms_relation(int n, int digits, double x) {
return -max(0.0, x * log10(M_E)) - (n + 1) * log10(abs(x))
+ (n + 1) * log10((n + 1) / M_E)
+ 0.5 * log10(2 * M_PI * (n + 1))
+ x * log10(M_E) - digits;
}
/**
* @brief Estimates the number of terms needed in a power series to reach a target precision.
*
* Uses binary search to find the minimum number of terms n such that the
* specified error function (typically based on Stirling’s approximation) exceeds
* the desired number of decimal digits at point x.
*
* @param digits The target number of decimal digits precision.
* @param x The evaluation point for the power series.
* @param func A function that estimates the error given (n, digits, x).
* Defaults to digits_to_terms_relation().
* @param low The lower bound of the search range (default 20).
* @param high The upper bound of the search range (default 1e9).
* @return The smallest n such that func(n, digits, x) > 0.
*
* @note The function assumes that func is monotonic in n
* and that it transitions from negative to positive.
*/
int estimate_terms(int digits, double x,
std::function<double(int, int, double)> func = digits_to_terms_relation,
int low = 20, int high = 1e9) {
while (high - low > 1) {
int mid = (low + high) / 2;
double y = func(mid, digits, x);
if (y > 0) {
high = mid;
} else {
low = mid;
}
}
return high;
}
/**
* @brief Compute \f(e^x\f).
*
* Function to compute factorial \f(e^x\f) using a Maclauren series.
*
* @param x Exp argument.
* @return double \f(e^x\f)
*/
mpreal maclauren_exp(const mpq_class& uv) {
int terms = estimate_terms(DIGITS, uv.get_d()) + 2;
cout << "Terms: " << terms << endl;
mpz_class P, Q;
binary_split(0, terms, uv, P, Q);
mpq_class result_uv{P, Q};
return result_uv.get_mpq_t();
}
/**
* @brief Main function.
*
* The entry point of the program.
*
* Syntax: maclauren_exp_bs [x] [num_digits]
* x can be a fraction like 2/3 or in decimal notation 1.5
*
* @return int Returns 0 on successful execution.
*/
int main(int argc, char **argv) {
// Get the number of DIGITS if present. Use strtod() so that we can use
// exponential notation like 1e6.
if (argc > 2) {
char* endptr;
DIGITS = static_cast<long>(strtod(argv[2], &endptr));
if (*endptr != '\0' || endptr == argv[2]) {
cerr << "Error: Invalid number format\n";
return 1;
}
}
// The value for x as a rational number u/v if present.
// We can use decimal notation like 0.1 or a fraction like 1/3.
mpq_class uv{"1"};
if (argc > 1) {
if (strchr(argv[1], '/') == nullptr) {
mpf_class x{argv[1]};
uv = x;
} else {
uv = argv[1];
}
}
// Precision in bits (~3.3 bits per decimal digit)
PRECISION = static_cast<long>(DIGITS * 3.32193 + 16);
mpreal::set_default_prec(PRECISION); // Set default precision in bits
// Use the high_resolution_clock to time execution
auto start = chrono::high_resolution_clock::now();
mpreal exp = maclauren_exp(uv);
auto stop = chrono::high_resolution_clock::now();
chrono::duration<double> elapsed = stop - start;
cout << "x = u/v = " << uv << endl;
cout << "e^x = " << exp.toString(DIGITS) << endl;
cout << "std::exp(x) = " << setprecision(16) << std::exp(uv.get_d()) << endl;
cout << "Computation Time: " << setprecision(3) << elapsed.count() << " seconds" << endl;
return 0;
}
x = 1
Terms: 13524
e^x = 2.718281828459045235360287471352662497757247093699959574966967627724076630
35354759457138217852516642742746639193200305992181741359662904357290033429526059
56307381323286279434907632338298807531952510190115738341879307021540891499348841
...
14492534279922013463805383423696437674288625951401461782018107341005654667082368
54312816339049676558789901487477972479202502227218169405159042170892104287552188
65830860845270842392865259753614629003778016700165467168160534329290757303146656
2485802
std::exp(x) = 2.718281828459045
Computation Time: 26.109585792 seconds
We can further enhance the performance of the binary splitting version of our program by optimizing the factorial calculation using memoization. Memoization is an optimization technique that accelerates computations by caching the results of expensive calls to pure functions. When the same input is encountered again, the previously computed result is returned from the cache instead of recomputing it.
Earlier, we found that mpreal::pow() is just as efficient as our custom ipow() function, which uses exponentiation by squaring. As a result, there’s little room for further optimization in how we compute the numerator of each term.
maclauren_exp_bs_opt.cpp
/**
* @file maclauren_exp_bs_opt.cpp
* @brief Program to compute \f(e^x\f).
*
* This program can compute \f(e^x\f) using the Maclauren series with the
* GNU Multiple Precision (GMP) library.
* We will use binary splitting algorithm to potentially
* speed up the computation.
*
* Syntax: maclauren_exp_bs_opt [x] [num_digits]
* x can be a fraction like 2/3 or in decimal notation 1.5
*
* get mpreal.h at https://github.com/advanpix/mpreal.git
*
* Compile with command:\n
* macOS:\n
* install homebrew https://brew.sh\n
* put mpreal.h in /opt/homebrew/include\n
* % brew install gcc\n
* % brew install gmp\n
* % brew install mpfr\n
* % g++ -std=c++20 -O3 -I/opt/homebrew/include -L/opt/homebrew/lib -o maclauren_exp_bs_opt maclauren_exp_bs_opt.cpp -lmpfr -lgmp
*
* Linux:\n
* put mpreal.h in /urs/local/include\n
* % sudo apt install libgmp-dev\n
* % sudo apt install libmpfr-dev\n
* % g++ -std=c++20 -O3 -o maclauren_exp_bs_opt maclauren_exp_bs_opt.cpp -lmpfr -lgmp
* @author Richard Lesh
* @date 2025-07-27
* @version 1.0
*/
#include <chrono>
#include <cmath>
#include <gmpxx.h> // GMP C++ Wrapper
#include <iomanip>
#include <iostream>
#include <mpreal.h> // MPFR C++ Wrapper
using namespace mpfr;
using namespace std;
long DIGITS = 16; // in decimal digits
long PRECISION; // in bits
using namespace std;
/**
* @brief Computes the factorial of a non-negative integer.
*
* This function recursively computes the factorial \f(n!\f) of a given
* non-negative integer `n`. The factorial of 0 and 1 is defined as 1. For
* values greater than 1, the function returns `n * factorial(n - 1)`.
* Speeds up computation by memoizing previous results.
*
* @param n A non-negative integer whose factorial is to be computed.
* @return The factorial of the input number `n`. Returns 1 if `n` is 0 or 1.
*
* @note The function does not overflow because arbitrary-precision integers are used.
*/
// NOLINTNEXTLINE(misc-no-recursion)
mpz_class factorial(long n) {
static vector<mpz_class> factorials_cache{1,1};
if (n <= 1) return 1;
if (factorials_cache.size() < n + 1) {
factorials_cache.emplace_back( n * factorial(n - 1));
}
return factorials_cache[n];
}
/**
* @brief Compute fast exponentiation.
*
* Function to compute \f(x^n\f) quickly using Exponentiation by Squaring.
*
* @param x Base argument
* @param n Power argument
* @return long \f(x^n\f)
*/
mpz_class ipow(const mpz_class& x, int n) {
// Use exponentiation by squaring
mpz_class x2 = x;
mpz_class result{1};
while (n > 0) {
if (n % 2 == 1) result *= x2;
x2 *= x2;
n /= 2;
}
return result;
}
/**
* @brief Recursive binary splitting algorithm for computing a rational sum.
*
* This function recursively computes a pair of integers (P, Q) such that
* \f[
* \frac{P}{Q} = \sum_{k=a}^{b-1} \frac{x^k}{k!}
* \f]
* using a binary splitting strategy for improved performance and numerical stability.
*
* @param a Lower bound of the summation interval (inclusive).
* @param b Upper bound of the summation interval (exclusive).
* @param x Rational base of the power series (as an `mpq_class`).
* @param[out] P Resulting numerator of the partial sum, accumulated via binary splitting.
* @param[out] Q Resulting denominator of the partial sum, accumulated via binary splitting.
*
* @note The result \f$ \frac{P}{Q} \f$ represents the sum of the terms
* \f$ \frac{x^k}{k!} \f$ for \f$ k \in [a, b) \f$.
*
* @warning Assumes that `factorial(a)` is computable for the given range.
* May not perform well for very large ranges without optimizations.
*/
// NOLINTNEXTLINE(misc-no-recursion)
void binary_split(int a, int b, const mpq_class& x, mpz_class& P, mpz_class& Q) {
if (b - a == 1) {
P = ipow(x.get_num(), a);
Q = ipow(x.get_den(), a) * factorial(a);
} else {
int m = (a + b) / 2;
mpz_class P1, P2;
mpz_class Q1, Q2;
binary_split(a, m, x, P1, Q1);
binary_split(m, b, x, P2, Q2);
P = P1 * Q2 + P2 * Q1;
Q = Q1 * Q2;
}
}
/**
* @brief Approximates the number of digits contributed by a term in a power series expansion.
*
* This function models the relationship between the number of terms n,
* desired decimal digits of precision, and the evaluation point x
* using Stirling’s approximation. It returns the difference between the estimated
* number of digits the partial sum contributes and the required number of digits.
*
* The approximation is based on:
* \f[
* f(n, D, x) = -max(0,x \cdot log_{10}(e)) - (n+1)\log_{10}(\left| x \right|)
* + (n + 1) ( \log_{10}(\frac{n + 1}{e}) + \frac{1}{2} \log_{10}(2\pi (n + 1))
* + x\log_{10}(e) - D
* \f]
*
* @param n The number of terms used in the approximation.
* @param digits The target number of decimal digits precision.
* @param x The evaluation point of the series.
* @return Estimated log10 error contribution; positive means more digits
* than needed.
*
* @note This function is intended to be used as a scoring or objective function
* in root-finding or binary search for optimal \p n.
*/
double digits_to_terms_relation(int n, int digits, double x) {
return -max(0.0, x * log10(M_E)) - (n + 1) * log10(abs(x))
+ (n + 1) * log10((n + 1) / M_E)
+ 0.5 * log10(2 * M_PI * (n + 1))
+ x * log10(M_E) - digits;
}
/**
* @brief Estimates the number of terms needed in a power series to reach the
* target precision.
*
* Uses binary search to find the minimum number of terms n such that the
* specified error function (typically based on Stirling’s approximation) exceeds
* the desired number of decimal digits at point x.
*
* @param digits The target number of decimal digits precision.
* @param x The evaluation point for the power series.
* @param func A function that estimates the error given (n, digits, x).
* Defaults to digits_to_terms_relation().
* @param low The lower bound of the search range (default 20).
* @param high The upper bound of the search range (default 1e9).
* @return The smallest n such that func(n, digits, x) > 0.
*
* @note The function assumes that func is monotonic in n
* and that it transitions from negative to positive.
*/
int estimate_terms(int digits, double x,
std::function<double(int, int, double)> func = digits_to_terms_relation,
int low = 20, int high = 1e9) {
while (high - low > 1) {
int mid = (low + high) / 2;
double y = func(mid, digits, x);
if (y > 0) {
high = mid;
} else {
low = mid;
}
}
return high;
}
/**
* @brief Compute \f(e^x\f).
*
* Function to compute factorial \f(e^x\f) using a Maclauren series.
*
* @param x Exp argument.
* @return double \f(e^x\f)
*/
mpreal maclauren_exp(const mpq_class& uv) {
int terms = estimate_terms(DIGITS, uv.get_d()) + 2;
cout << "Terms: " << terms << endl;
mpz_class P, Q;
binary_split(0, terms, uv, P, Q);
mpq_class result_uv{P, Q};
return result_uv.get_mpq_t();
}
/**
* @brief Main function.
*
* The entry point of the program.
*
* Syntax: maclauren_exp_bs [x] [num_digits]
* x can be a fraction like 2/3 or in decimal notation 1.5
*
* @return int Returns 0 on successful execution.
*/
int main(int argc, char **argv) {
// Get the number of DIGITS if present. Use strtod() so that we can use
// exponential notation like 1e6.
if (argc > 2) {
char* endptr;
DIGITS = static_cast<long>(strtod(argv[2], &endptr));
if (*endptr != '\0' || endptr == argv[2]) {
cerr << "Error: Invalid number format\n";
return 1;
}
}
// The value for x as a rational number u/v if present.
// We can use decimal notation like 0.1 or a fraction like 1/3.
mpq_class uv{"1"};
if (argc > 1) {
if (strchr(argv[1], '/') == nullptr) {
mpf_class x{argv[1]};
uv = x;
} else {
uv = argv[1];
}
}
// Precision in bits (~3.3 bits per decimal digit)
PRECISION = static_cast<long>(DIGITS * 3.32193 + 16);
mpreal::set_default_prec(PRECISION); // Set default precision in bits
// Use the high_resolution_clock to time execution
auto start = chrono::high_resolution_clock::now();
mpreal exp = maclauren_exp(uv);
auto stop = chrono::high_resolution_clock::now();
chrono::duration<double> elapsed = stop - start;
cout << "x = u/v = " << uv << endl;
cout << "e^x = " << exp.toString(DIGITS) << endl;
cout << "std::exp(x) = " << setprecision(16) << std::exp(uv.get_d()) << endl;
cout << "Computation Time: " << setprecision(3) << elapsed.count() << " seconds" << endl;
return 0;
}
When comparing the performance of our binary splitting implementations, we observe that the binary split version performs significantly worse than the unoptimized direct Maclaurin computation. Even with optimizations, the binary split approach still lags behind the optimized Maclaurin version. This is due to the additional overhead introduced by the binary splitting algorithm.
Binary splitting only becomes advantageous when computing a large number of significant digits, where its efficiency begins to outweigh the setup cost. Additionally, the terms in the \(e^x\) series are relatively simple to compute, offering little benefit from the binary split method. This technique is far more effective for series with more computationally intensive terms, such as the Chudnovsky series used for calculating π.
Additionally, because binary splitting defers division operations, it can generate extremely large numerators and denominators. This can result in a great deal of memory pressure.
Advanced Binary Splitting
None of the algorithms discussed so far are efficient enough to compute a million or more digits within a reasonable timeframe. A more advanced method, based on binary splitting, was introduced by Bruno Haible and Thomas Papanikolaou in their 1998 paper Fast Multiprecision Evaluation of Series of Rational Numbers. Similar to our optimized Maclaurin implementation, this approach leverages the fact that each term in the series can be efficiently derived from previous terms. When combined with the binary splitting technique and some multi-threading, it yields an exceptionally fast algorithm.
binary_split_exp.cpp
/**
* @file binary_split_e.cpp
* @brief Program to compute \(e^x\).
*
* This program can compute \(e^x\) using the binary
* splitting technique of https://www.ginac.de/CLN/binsplit.pdf
* Requires ThreadPool class.
*
* Syntax: binary_split_e [x] [num_digits]
*
* get mpreal.h at https://github.com/advanpix/mpreal.git
*
* Compile with command:\n
* macOS:\n
* install homebrew https://brew.sh\n
* put mpreal.h in /opt/homebrew/include\n
* brew install gcc\n
* brew install gmp\n
* brew install mpfr\n
* g++ -std=c++20 -O3 -I/opt/homebrew/include -L/opt/homebrew/lib -o binary_split_exp binary_split_exp.cpp ThreadPool.cpp -lmpfr -lgmp
*
* Linux\n
* put mpreal.h in /urs/local/include\n
* sudo apt install libgmp-dev\n
* sudo apt install libmpfr-dev\n
* g++ -std=c++20 -O3 -o binary_split_exp binary_split_exp.cpp ThreadPool.cpp -lmpfr -lgmp
*
* @author Richard Lesh
* @date 2025-07-26
* @version 1.0
*/
#include <cassert>
#include <cmath>
#include <gmpxx.h> // GMP C++ Wrapper
#include <iostream>
#include <iomanip>
#include <mpreal.h> // MPFR C++ Wrapper
#include <ThreadPool.hpp>
#include <semaphore>
using namespace mpfr;
using namespace std;
long DIGITS = 50; // in decimal digits
long PRECISION; // in bits
long NUM_TERMS;
int MAX_THREAD_DEPTH = static_cast<int>(round(log(thread::hardware_concurrency()) / log(2.0)));
struct PQBT {
mpz_class P;
mpz_class Q;
mpz_class B;
mpz_class T;
};
/**
* @brief Compute e.
*
* Function to compute \f( e^{\frac{u}{v}} \f) using binary splitting.
*
* @return mpreal A multi-precision value for e.
*/
// NOLINTNEXTLINE(misc-no-recursion)
PQBT binary_split(long n, long m, int depth, const mpz_class& mpz_u, const mpz_class& mpz_v) {
PQBT pqbt;
if (m - n == 1) {
if (n == 0) {
pqbt.P = 1; // p(0)
pqbt.Q = 1; // q(0)
pqbt.B = 1; // b(0)
pqbt.T = 1; // a(0)p(0)
} else {
// p(n)
pqbt.P = mpz_u;
// q(n)
pqbt.Q = mpz_v * n;
// b(n)
pqbt.B = 1;
// a(n)p(n)
pqbt.T = mpz_u ;
}
return pqbt;
} else if (m - n == 2) {
PQBT left = binary_split(n, m - 1, depth + 1, mpz_u, mpz_v);
PQBT right = binary_split(n + 1, m, depth + 1, mpz_u, mpz_v);
pqbt.P = left.P * right.P;
pqbt.Q = left.Q * right.Q;
pqbt.B = left.B * right.B;
pqbt.T = right.B * right.Q * left.T + left.B * left.P * right.T;
return pqbt;
} else {
long new_m = (m + n) / 2;
PQBT left, right;
if (depth < MAX_THREAD_DEPTH) {
shared_ptr<ThreadPool> pool = ThreadPool::getThreadPool();
counting_semaphore<1> semaphore(0);
pool->enqueue([&right, &pool, &semaphore, &mpz_u, &mpz_v,
m, new_m, depth]() {
pool->setThreadName(to_string(new_m) + "-" + to_string(m));
right = binary_split(new_m, m, depth + 1, mpz_u, mpz_v);
semaphore.release();
});
left = binary_split(n, new_m, depth + 1, mpz_u, mpz_v);
semaphore.acquire();
} else {
right = binary_split(new_m, m, depth + 1, mpz_u, mpz_v);
left = binary_split(n, new_m, depth + 1, mpz_u, mpz_v);
}
pqbt.P = left.P * right.P;
pqbt.Q = left.Q * right.Q;
pqbt.B = left.B * right.B;
pqbt.T = right.B * right.Q * left.T + left.B * left.P * right.T;
return pqbt;
}
}
// Function to solve
double digits_function(int n, int digits, double x) {
return -max(0.0, x * log10(M_E)) - (n + 1) * log10(abs(x))
+ (n + 1) * log10((n + 1) / M_E)
+ 0.5 * log10(2 * M_PI * (n + 1))
+ x * log10(M_E) - digits;
}
int estimate_terms(int digits, double x,
std::function<double(int, int, double)> func = digits_function,
int low = 1, int high = 1e9) {
while (high - low > 1) {
int mid = (low + high) / 2;
double y = func(mid, digits, x);
if (y > 0) {
high = mid;
} else {
low = mid;
}
}
return high;
}
mpreal exp(const mpq_class& x) {
mpq_class uv{x};
NUM_TERMS = estimate_terms(DIGITS, uv.get_d()) + 2;
mpz_class mpz_u{1}, mpz_v{1};
PQBT sum = binary_split(0, NUM_TERMS, 0, uv.get_num(), uv.get_den());
mpz_class d = sum.B * sum.Q;
mpreal::set_default_prec(PRECISION); // Set default precision in bits
long num_digits_t = mpz_sizeinbase(sum.T.get_mpz_t(), 2);
long num_digits_d = mpz_sizeinbase(d.get_mpz_t(), 2);
mpfr::mpreal::set_emin(- std::max(num_digits_t, num_digits_d) - 2);
mpfr::mpreal::set_emax(std::max(num_digits_t, num_digits_d) + 2);
mpreal num{sum.T.get_mpz_t()};
mpreal denom{d.get_mpz_t()};
mpreal exp = num / denom;
return exp;
}
/**
* @brief Main function.
*
* The entry point of the program.
*
* Syntax: binary_split_exp [x] [num_digits]
*
* @return int Returns 0 on successful execution.
*/
int main(int argc, char **argv) {
// Get the number of DIGITS if present. Use strtod() so that we can use
// exponential notation like 1e6.
if (argc > 2) {
char* endptr;
DIGITS = static_cast<long>(strtod(argv[2], &endptr));
if (*endptr != '\0' || endptr == argv[2]) {
cerr << "Error: Invalid number format\n";
return 1;
}
}
// The value for x as a rational number u/v if present.
// We can use decimal notation like 0.1 or a fraction like 1/3.
mpq_class uv{"1"};
if (argc > 1) {
if (strchr(argv[1], '/') == nullptr) {
mpf_class x{argv[1]};
uv = x;
} else {
uv = argv[1];
}
}
// Precision in bits (~3.3 bits per decimal digit)
PRECISION = static_cast<long>(DIGITS * 3.32193 + 16);
auto start = chrono::high_resolution_clock::now();
mpreal e = exp(uv);
auto stop = chrono::high_resolution_clock::now();
assert(!isnan(e));
chrono::duration<double> elapsed = stop - start;
cout << "x = u/v = " << uv << endl;
cout << "e^x = " << e.toString(DIGITS) << endl;
cout << "std::exp(x) = " << setprecision(16) << std::exp(uv.get_d()) << endl;
cout << "Computation Time: " << setprecision(3) << elapsed.count() << " seconds" << endl;
cout << "Terms: " << NUM_TERMS << endl;
return 0;
}
Time to Compute \(e^1\)
| Digits | MP | MP Opt | BS | BS Opt | Adv. Binary Split |
|---|---|---|---|---|---|
| 1E+02 | 0.000155 | 4.1E-05 | 0.000115 | 7.46E-05 | 0.000521 |
| 5E+02 | 0.00175 | 0.000341 | 0.0023 | 0.0016 | 0.000390 |
| 1E+03 | 0.00683 | 0.00145 | 0.0112 | 0.0082 | 0.000389 |
| 5E+03 | 0.204 | 0.0585 | 0.458 | 0.354 | 0.000469 |
| 1E+04 | 0.935 | 0.387 | 2.15 | 1.67 | 0.000531 |
| 2E+04 | 4.31 | 1.61 | 10.2 | 8.06 | 0.000827 |
| 3E+04 | 11.3 | 4.25 | 25.8 | 20 | 0.001104 |
| 4E+04 | 23.3 | 8.17 | 51.1 | 37.8 | 0.001503 |
| 5E+04 | 42 | 14.2 | 86.2 | 61.3 | 0.001805 |
| 6E+04 | 62.9 | 21 | 131 | 91.3 | 0.002223 |
| 7E+04 | 95.1 | 28.7 | 187 | 126 | 0.002444 |
| 8E+04 | 136 | 40.2 | 260 | 172 | 0.002716 |
| 9E+04 | 181 | 51.9 | 337 | 219 | 0.003162 |
| 1E+05 | 237 | 65.4 | 428 | 267 | 0.003675 |
| 1.25E+05 | 408 | 104 | 743 | 447 | 0.005059 |
| 1.5E+05 | 699 | 165 | 1160 | 660 | 0.005938 |
| 1.75E+05 | 1040 | 229 | 1690 | 912 | 0.007357 |
| 2E+05 | 1440 | 298 | 2310 | 1190 | 0.007484 |
| 3E+05 | 4240 | 700 | 6390 | 2980 | 0.011925 |
| 5E+05 | 16500 | 1900 | 23700 | 8820 | 0.019893 |
We can see that the simple binary split programs (BS) perform worse than the direct series computation versions (MP). This is because the binary splitting introduces overhead that must be overcome. You can see in the graphs below that the optimized BS program eventually out performs the naive direct series algorithm. Even so, all four of our first four programs have very similar computational complexity close to \(O(n^2)\). It isn't until we combine all these techniques in the advanced binary split program that we see a change in the computation complexity closer to \(O(n \cdot \log(n))\).
Reduction of exp(x)
One final optimization that we can apply for large values of \(x\), is to reduce the problem to a simpler problem. Since the series for \(e^x\) converges faster near \(x = 0\), we can improve convergence by transforming the computation of \(e^x\) into the computation of \(e^y\) for some \(y \leq 1\). This is done by repeatedly dividing \(x\) by 2 until \(y = \frac{x}{2^n} \leq 1\), where n is the number of halvings. Then, noting that \(x = 2^n \cdot y\), we have:
\[ e^x = e^{2^n \cdot y} = (e^y)^{2^n}\]
We also note that:
\[ e^{-x} = \frac{1}{e^x} \]
So we can also transform computations with large negative \(x\) into simpler computations.
binary_split_exp_opt.cpp
/**
* @file binary_split_exp_opt.cpp
* @brief Program to compute \(e^x\).
*
* This program can compute \(e^x\) using the binary
* splitting technique of https://www.ginac.de/CLN/binsplit.pdf
* Requires ThreadPool class.
*
* Syntax: binary_split_e [x] [num_digits]
*
* get mpreal.h at https://github.com/advanpix/mpreal.git
*
* Compile with command:\n
* macOS:\n
* install homebrew https://brew.sh\n
* put mpreal.h in /opt/homebrew/include\n
* brew install gcc\n
* brew install gmp\n
* brew install mpfr\n
* g++ -std=c++20 -O3 -I/opt/homebrew/include -L/opt/homebrew/lib -o binary_split_exp_opt binary_split_exp_opt.cpp ThreadPool.cpp -lmpfr -lgmp
*
* Linux\n
* put mpreal.h in /urs/local/include\n
* sudo apt install libgmp-dev\n
* sudo apt install libmpfr-dev\n
* g++ -std=c++20 -O3 -o v binary_split_exp_opt.cpp ThreadPool.cpp -lmpfr -lgmp
*
* @author Richard Lesh
* @date 2025-07-26
* @version 1.0
*/
#include <cassert>
#include <cmath>
#include <gmpxx.h> // GMP C++ Wrapper
#include <iostream>
#include <iomanip>
#include <mpreal.h> // MPFR C++ Wrapper
#include <ThreadPool.hpp>
#include <semaphore>
using namespace mpfr;
using namespace std;
long DIGITS = 50; // in decimal digits
long PRECISION; // in bits
long NUM_TERMS;
int MAX_THREAD_DEPTH = static_cast<int>(round(log(thread::hardware_concurrency()) / log(2.0)));
struct PQBT {
mpz_class P;
mpz_class Q;
mpz_class B;
mpz_class T;
};
/**
* @brief Compute e.
*
* Function to compute \f( e^{\frac{u}{v}} \f) using binary splitting.
*
* @return mpreal A multi-precision value for e.
*/
// NOLINTNEXTLINE(misc-no-recursion)
PQBT binary_split(long n, long m, int depth, const mpz_class& mpz_u, const mpz_class& mpz_v) {
PQBT pqbt;
if (m - n == 1) {
if (n == 0) {
pqbt.P = 1; // p(0)
pqbt.Q = 1; // q(0)
pqbt.B = 1; // b(0)
pqbt.T = 1; // a(0)p(0)
} else {
// p(n)
pqbt.P = mpz_u;
// q(n)
pqbt.Q = mpz_v * n;
// b(n)
pqbt.B = 1;
// a(n)p(n)
pqbt.T = mpz_u ;
}
return pqbt;
} else if (m - n == 2) {
PQBT left = binary_split(n, m - 1, depth + 1, mpz_u, mpz_v);
PQBT right = binary_split(n + 1, m, depth + 1, mpz_u, mpz_v);
pqbt.P = left.P * right.P;
pqbt.Q = left.Q * right.Q;
pqbt.B = left.B * right.B;
pqbt.T = right.B * right.Q * left.T + left.B * left.P * right.T;
return pqbt;
} else {
long new_m = (m + n) / 2;
PQBT left, right;
if (depth < MAX_THREAD_DEPTH) {
shared_ptr<ThreadPool> pool = ThreadPool::getThreadPool();
counting_semaphore<1> semaphore(0);
pool->enqueue([&right, &pool, &semaphore, &mpz_u, &mpz_v,
m, new_m, depth]() {
pool->setThreadName(to_string(new_m) + "-" + to_string(m));
right = binary_split(new_m, m, depth + 1, mpz_u, mpz_v);
semaphore.release();
});
left = binary_split(n, new_m, depth + 1, mpz_u, mpz_v);
semaphore.acquire();
} else {
right = binary_split(new_m, m, depth + 1, mpz_u, mpz_v);
left = binary_split(n, new_m, depth + 1, mpz_u, mpz_v);
}
pqbt.P = left.P * right.P;
pqbt.Q = left.Q * right.Q;
pqbt.B = left.B * right.B;
pqbt.T = right.B * right.Q * left.T + left.B * left.P * right.T;
return pqbt;
}
}
// Function to solve
double digits_function(int n, int digits, double x) {
return -max(0.0, x * log10(M_E)) - (n + 1) * log10(abs(x))
+ (n + 1) * log10((n + 1) / M_E)
+ 0.5 * log10(2 * M_PI * (n + 1))
+ x * log10(M_E) - digits;
}
int estimate_terms(int digits, double x,
std::function<double(int, int, double)> func = digits_function,
int low = 1, int high = 1e9) {
while (high - low > 1) {
int mid = (low + high) / 2;
double y = func(mid, digits, x);
if (y > 0) {
high = mid;
} else {
low = mid;
}
}
return high;
}
mpreal exp(const mpq_class& x) {
mpq_class uv{x};
// Reduce problems with input -x to 1/x
if (sgn(uv) < 0) {
mpq_class uv_neg{abs(uv)};
return mpreal(1.0) / exp(uv_neg);
}
// Reduce the problem to one where x is closer to 0.
long factors = 1;
while (uv.get_num() > uv.get_den()) {
uv = uv / 2;
factors *= 2;
}
NUM_TERMS = estimate_terms(DIGITS, uv.get_d()) + 2;
mpz_class mpz_u{1}, mpz_v{1};
PQBT sum = binary_split(0, NUM_TERMS, 0, uv.get_num(), uv.get_den());
mpz_class d = sum.B * sum.Q;
mpreal::set_default_prec(PRECISION); // Set default precision in bits
long num_digits_t = mpz_sizeinbase(sum.T.get_mpz_t(), 2);
long num_digits_d = mpz_sizeinbase(d.get_mpz_t(), 2);
mpfr::mpreal::set_emin(- std::max(num_digits_t, num_digits_d) - 2);
mpfr::mpreal::set_emax(std::max(num_digits_t, num_digits_d) + 2);
mpreal num{sum.T.get_mpz_t()};
mpreal denom{d.get_mpz_t()};
mpreal exp = num / denom;
if (factors > 1) return pow(exp, factors);
return exp;
}
/**
* @brief Main function.
*
* The entry point of the program.
*
* Syntax: binary_split_exp_opt [x] [num_digits]
*
* @return int Returns 0 on successful execution.
*/
int main(int argc, char **argv) {
// Get the number of DIGITS if present. Use strtod() so that we can use
// exponential notation like 1e6.
if (argc > 2) {
char* endptr;
DIGITS = static_cast<long>(strtod(argv[2], &endptr));
if (*endptr != '\0' || endptr == argv[2]) {
cerr << "Error: Invalid number format\n";
return 1;
}
}
// The value for x as a rational number u/v if present.
// We can use decimal notation like 0.1 or a fraction like 1/3.
mpq_class uv{"1"};
if (argc > 1) {
if (strchr(argv[1], '/') == nullptr) {
mpf_class x{argv[1]};
uv = x;
} else {
uv = argv[1];
}
}
// Precision in bits (~3.3 bits per decimal digit)
PRECISION = static_cast<long>(DIGITS * 3.32193 + 16);
auto start = chrono::high_resolution_clock::now();
mpreal e = exp(uv);
auto stop = chrono::high_resolution_clock::now();
assert(!isnan(e));
chrono::duration<double> elapsed = stop - start;
cout << "x = u/v = " << uv << endl;
cout << "e^x = " << e.toString(DIGITS) << endl;
cout << "std::exp(x) = " << setprecision(16) << std::exp(uv.get_d()) << endl;
cout << "Computation Time: " << setprecision(3) << elapsed.count() << " seconds" << endl;
cout << "Terms: " << NUM_TERMS << endl;
return 0;
}
After applying this optimization we find that for large values of \(x\) it significantly reduces the number of terms needed to compute \(e^x\).
Computing π
Using the Haible and Papanikolaou paper we also find a very fast algorithm to compute π using binary splitting and the very fast but complex Chudnovsky series for π. This is a good candidate for binary splitting because each term in the series it costly to compute. See the post Computing π (Revisited) for details on how to apply binary splitting to the computation of π.
Computing \(\sqrt{x}\)
To compute the square root of x we need a different technique. It is an iterative root-finding algorithm called Newton–Raphson. It produces successively better approximations to the roots (or zeroes) of a real-valued function. All it requires is a function \(f\), its derivitive \(f’\) and an initial guess \(x_0\). To compute the next better approximation we use:
\[x_1 = x_0 - \frac{f(x_0)}{f'(x_0)}\]
Geometrically, \((x_1, 0)\) is the x-intercept of the tangent to \(f\) at \((x_0, f(x_0)\), e.g., the improved guess, x1, is the unique root of the linear approximation of \(f\) at the initial guess, \(x_0\) In general we have
\[x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}\]
until a sufficiently precise value is reached. When the iterations get close to the solution, the number of correct digits roughly doubles with each step.
So to compute the square root of N we use the function \(f(x) = x^2 - N\) because \(\sqrt{N}\) is the root, e.g. solution to \(f(x) = 0\). The derivative of this function is \(f’(x) = 2x\). So this gives us the iterative formula:
\[x_{n+1} = x_n - \frac{x_n^2 - N}{2x_n} = \frac{x_n^2 + N}{2x_n}\]
square_root.cpp
/**
* @file square_root.cpp
* @brief Program to compute square roots.
*
* This program can compute the square root of N using
* the Newton–Raphson method.
*
* Syntax: square_root [N] [num_digits]
*
* get mpreal.h at https://github.com/advanpix/mpreal.git
*
* Compile with command:\n
* macOS:\n
* install homebrew https://brew.sh\n
* put mpreal.h in /opt/homebrew/include\n
* brew install gcc\n
* brew install gmp\n
* brew install mpfr\n
* g++ -std=c++20 -O3 -I/opt/homebrew/include -L/opt/homebrew/lib -o square_root square_root.cpp -lmpfr -lgmp
*
* Linux\n
* put mpreal.h in /urs/local/include\n
* sudo apt install libgmp-dev\n
* sudo apt install libmpfr-dev\n
* g++ -std=c++20 -O3 -o square_root square_root.cpp -lmpfr -lgmp
*
* @author Richard Lesh
* @date 2025-07-26
* @version 1.0
*/
#include <cassert>
#include <cmath>
#include <gmpxx.h> // GMP C++ Wrapper
#include <iostream>
#include <iomanip>
#include <mpreal.h> // MPFR C++ Wrapper
using namespace mpfr;
using namespace std;
long DIGITS = 50; // in decimal digits
long PRECISION; // in bits
/**
* @brief Compute square root.
*
* Function to compute square root using Newton–Raphson.
*
* @param N Square root argument.
* @param num_digits Accuracy in decimal digits
* @return mpreal A multi-precision value for \f( \sqrt{N} \f).
*/
mpreal sqrt(const mpreal &N, long num_digits) {
if (N < 0) throw domain_error("Square root undefined for N < 0");
if (N == 0) return 0;
mpreal x = N / 4.0;
mpreal max_relative_err = pow(10, -num_digits);
int i = 0;
mpreal last = x;
do {
last = x;
x = (x * x + N) / (2 * x);
++i;
} while (abs(x - last) / x > max_relative_err);
cout << "Iterations: " << i << endl;
return x;
}
/**
* @brief Main function.
*
* The entry point of the program.
*
* Syntax: square_root [N] [num_digits]
*
* @return int Returns 0 on successful execution.
*/
int main(int argc, char **argv) {
if (argc > 2) {
char* endptr;
DIGITS = static_cast<long>(strtod(argv[2], &endptr));
if (*endptr != '\0' || endptr == argv[2]) {
cerr << "Error: Invalid number format\n";
return 1;
}
}
// Precision in bits (~3.3 bits per decimal digit)
PRECISION = static_cast<long>(DIGITS * 3.32193 + 16);
mpreal::set_default_prec(PRECISION); // Set default precision in bits
mpfr::mpreal::set_emin(-PRECISION - 2);
mpfr::mpreal::set_emax(PRECISION + 2);
mpreal x = 2.0;
if (argc > 1) {
x = mpreal(argv[1]);
}
auto start = chrono::high_resolution_clock::now();
mpreal sqrt = ::sqrt(x, DIGITS);
auto stop = chrono::high_resolution_clock::now();
assert(!isnan(sqrt));
chrono::duration<double> elapsed = stop - start;
cout << "x = " << x << endl;
cout << "sqrt(x) = " << sqrt.toString(DIGITS) << endl;
cout << "std::sqrt(x) = " << setprecision(16) << std::sqrt(x.toDouble()) << endl;
printf("Computation Time: %f seconds\n", elapsed.count());
return 0;
}
Computing \(\sqrt[3]{x}\) and higher
We use the same approach to compute the cube root as we did to compute the square root. We use Newton-Pahpson method to find the root of the equation:
\[f(x) = x^3 - N = 0\]
Where \(N\) is the value we want to use to find the cube root. We use the function and its derivative to compute the iteration formula:
\[ x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)} = x_n - \frac{x_n^3 - N}{3x_n^2} = \frac{2x_n + \frac{N}{x_n^2}}{3} \]
cube_root.cpp
/**
* @file cube_root.cpp
* @brief Program to compute cube roots.
*
* This program can compute the cube root of N using
* the Newton–Raphson method.
*
* Syntax: cube_root [N] [num_digits]
*
* get mpreal.h at https://github.com/advanpix/mpreal.git
*
* Compile with command:\n
* macOS:\n
* install homebrew https://brew.sh\n
* put mpreal.h in /opt/homebrew/include\n
* brew install gcc\n
* brew install gmp\n
* brew install mpfr\n
* g++ -std=c++20 -O3 -I/opt/homebrew/include -L/opt/homebrew/lib -o cube_root cube_root.cpp -lmpfr -lgmp
*
* Linux\n
* put mpreal.h in /urs/local/include\n
* sudo apt install libgmp-dev\n
* sudo apt install libmpfr-dev\n
* g++ -std=c++20 -O3 -o cube_root cube_root.cpp -lmpfr -lgmp
*
* @author Richard Lesh
* @date 2025-07-26
* @version 1.0
*/
#include <cassert>
#include <cmath>
#include <gmpxx.h> // GMP C++ Wrapper
#include <iostream>
#include <iomanip>
#include <mpreal.h> // MPFR C++ Wrapper
using namespace mpfr;
using namespace std;
long DIGITS = 50; // in decimal digits
long PRECISION; // in bits
/**
* @brief Compute cube root.
*
* Function to compute cube root using Newton–Raphson.
*
* @param N Cube root argument.
* @param num_digits Accuracy in decimal digits
* @return mpreal A multi-precision value for \f( \sqrt[3]{N} \f).
*/
mpreal cbrt(const mpreal &N, long num_digits) {
if (N == 0) throw domain_error("Cube root undefined for N = 0");
mpreal x = N / 4.0;
mpreal max_relative_err = pow(10, -num_digits);
int i = 0;
mpreal last = x;
do {
last = x;
x = (2 * x + N / (x * x)) / 3;
++i;
} while (abs(x - last) / abs(x) > max_relative_err);
cout << "Iterations: " << i << endl;
return x;
}
/**
* @brief Main function.
*
* The entry point of the program.
*
* Syntax: cube_root [N] [num_digits]
*
* @return int Returns 0 on successful execution.
*/
int main(int argc, char **argv) {
if (argc > 2) {
char* endptr;
DIGITS = static_cast<long>(strtod(argv[2], &endptr));
if (*endptr != '\0' || endptr == argv[2]) {
cerr << "Error: Invalid number format\n";
return 1;
}
}
// Precision in bits (~3.3 bits per decimal digit)
PRECISION = static_cast<long>(DIGITS * 3.32193 + 16);
mpreal::set_default_prec(PRECISION); // Set default precision in bits
mpfr::mpreal::set_emin(-PRECISION - 2);
mpfr::mpreal::set_emax(PRECISION + 2);
mpreal x = 2.0;
if (argc > 1) {
x = mpreal(argv[1]);
}
auto start = chrono::high_resolution_clock::now();
mpreal cbrt = ::cbrt(x, DIGITS);
auto stop = chrono::high_resolution_clock::now();
assert(!isnan(cbrt));
chrono::duration<double> elapsed = stop - start;
cout << "x = " << x << endl;
cout << "cbrt(x) = " << cbrt.toString(DIGITS) << endl;
if (x > 0)
cout << "std::pow(x, 1/3) = " << setprecision(16) << std::pow(x.toDouble(), 1./3.) << endl;
else
cout << "std::pow(x, 1/3) = " << setprecision(16) << -std::pow(-x.toDouble(), 1./3.) << endl;
printf("Computation Time: %f seconds\n", elapsed.count());
return 0;
}
For the general \(nth\) root of N, \(\sqrt[n]{N}\) we have:
\[ x_{k+1} = \frac{1}{n} \left[ (n - 1)x_k + \frac{N}{x_k^{n - 1}} \right] \]
nth_root
/**
* @file nth_root.cpp
* @brief Program to compute nth roots.
*
* This program can compute the nth root of N using
* the Newton–Raphson method.
*
* Syntax: nth_root [N] [nth_root] [num_digits]
*
* get mpreal.h at https://github.com/advanpix/mpreal.git
*
* Compile with command:\n
* macOS:\n
* install homebrew https://brew.sh\n
* put mpreal.h in /opt/homebrew/include\n
* brew install gcc\n
* brew install gmp\n
* brew install mpfr\n
* g++ -std=c++20 -O3 -I/opt/homebrew/include -L/opt/homebrew/lib -o nth_root nth_root.cpp -lmpfr -lgmp
*
* Linux\n
* put mpreal.h in /urs/local/include\n
* sudo apt install libgmp-dev\n
* sudo apt install libmpfr-dev\n
* g++ -std=c++20 -O3 -o nth_root nth_root.cpp -lmpfr -lgmp
*
* @author Richard Lesh
* @date 2025-07-26
* @version 1.0
*/
#include <cassert>
#include <cmath>
#include <gmpxx.h> // GMP C++ Wrapper
#include <iostream>
#include <iomanip>
#include <mpreal.h> // MPFR C++ Wrapper
using namespace mpfr;
using namespace std;
long DIGITS = 50; // in decimal digits
long PRECISION; // in bits
/**
* @brief Compute nth root.
*
* Function to compute nth root using Newton–Raphson.
*
* @param N nth root argument.
* @param nth_root The root to compute.
* @param num_digits Accuracy in decimal digits
* @return mpreal A multi-precision value for \f( \sqrt[n]{N} \f).
*/
mpreal nth_root(const mpreal &N, int nth_root, long num_digits) {
if (nth_root < 1) throw domain_error("nth root undefined for n < 1");
if (nth_root % 2 == 0) {
if (N == 0) return 0;
else if (N < 0) throw domain_error("nth root undefined for N < 0");
} else {
if (N == 0) throw domain_error("nth root undefined for N = 0");
}
mpreal x = N / nth_root;
mpreal max_relative_err = pow(10, -num_digits);
int i = 0;
mpreal last = x;
do {
last = x;
x = ((nth_root - 1) * x + N / pow(x, nth_root - 1)) / nth_root;
++i;
} while (abs(x - last) / abs(x) > max_relative_err);
cout << "Iterations: " << i << endl;
return x;
}
/**
* @brief Main function.
*
* The entry point of the program.
*
* Syntax: nth_root [N] [nth_root] [num_digits]
*
* @return int Returns 0 on successful execution.
*/
int main(int argc, char **argv) {
if (argc > 3) {
char* endptr;
DIGITS = static_cast<long>(strtod(argv[3], &endptr));
if (*endptr != '\0' || endptr == argv[3]) {
cerr << "Error: Invalid number format\n";
return 1;
}
}
int n = 2;
if (argc > 2) {
n = atoi(argv[2]);
}
// Precision in bits (~3.3 bits per decimal digit)
PRECISION = static_cast<long>(DIGITS * 3.32193 + 16);
mpreal::set_default_prec(PRECISION); // Set default precision in bits
mpfr::mpreal::set_emin(-PRECISION - 2);
mpfr::mpreal::set_emax(PRECISION + 2);
mpreal x = 2.0;
if (argc > 1) {
x = mpreal(argv[1]);
}
auto start = chrono::high_resolution_clock::now();
mpreal nth_root = ::nth_root(x, n, DIGITS);
auto stop = chrono::high_resolution_clock::now();
assert(!isnan(nth_root));
chrono::duration<double> elapsed = stop - start;
cout << "x = " << x << endl;
cout << "nth_root(x, " << n << ") = " << nth_root.toString(DIGITS) << endl;
if (x > 0)
cout << "std::pow(x, 1/" << n << ") = " << setprecision(16) << std::pow(x.toDouble(), 1./n) << endl;
else if (n % 2 == 1)
cout << "std::pow(x, 1/" << n << ") = " << setprecision(16) << -std::pow(-x.toDouble(), 1./n) << endl;
printf("Computation Time: %f seconds\n", elapsed.count());
return 0;
}
/**
* @file nth_root.cpp
* @brief Program to nth roots.
*
* This program can compute the nth root of N using
* the Newton–Raphson method.
*
* Syntax: nth_root [N] [num_digits]
*
* get mpreal.h at https://github.com/advanpix/mpreal.git
*
* Compile with command:\n
* macOS:\n
* install homebrew https://brew.sh\n
* put mpreal.h in /opt/homebrew/include\n
* brew install gcc\n
* brew install gmp\n
* brew install mpfr\n
* g++ -std=c++20 -O3 -I/opt/homebrew/include -L/opt/homebrew/lib -o nth_root nth_root.cpp -lmpfr -lgmp
*
* Linux\n
* put mpreal.h in /urs/local/include\n
* sudo apt install libgmp-dev\n
* sudo apt install libmpfr-dev\n
* g++ -std=c++20 -O3 -o nth_root nth_root.cpp -lmpfr -lgmp
*
* @author Richard Lesh
* @date 2025-07-26
* @version 1.0
*/
#include <cassert>
#include <cmath>
#include <gmpxx.h> // GMP C++ Wrapper
#include <iostream>
#include <iomanip>
#include <mpreal.h> // MPFR C++ Wrapper
using namespace mpfr;
using namespace std;
long DIGITS = 50; // in decimal digits
long PRECISION; // in bits
/**
* @brief Compute nth root.
*
* Function to compute nth root using Newton–Raphson.
*
* @param N nth root argument.
* @param nth_root The root to compute.
* @param num_digits Accuracy in decimal digits
* @return mpreal A multi-precision value for \f( \sqrt[n]{N} \f).
*/
mpreal nth_root(const mpreal &N, int nth_root, long num_digits) {
if (nth_root < 1) throw domain_error("nth root undefined for n < 1");
if (nth_root % 2 == 0) {
if (N == 0) return 0;
else if (N < 0) throw domain_error("nth root undefined for N < 0");
} else {
if (N == 0) throw domain_error("nth root undefined for N = 0");
}
mpreal x = N / nth_root;
mpreal max_relative_err = pow(10, -num_digits);
int i = 0;
mpreal last = x;
do {
last = x;
x = ((nth_root - 1) * x + N / pow(x, nth_root - 1)) / nth_root;
++i;
} while (abs(x - last) / abs(x) > max_relative_err);
cout << "Iterations: " << i << endl;
return x;
}
/**
* @brief Main function.
*
* The entry point of the program.
*
* Syntax: nth_root [N] [num_digits]
*
* @return int Returns 0 on successful execution.
*/
int main(int argc, char **argv) {
if (argc > 3) {
char* endptr;
DIGITS = static_cast<long>(strtod(argv[3], &endptr));
if (*endptr != '\0' || endptr == argv[3]) {
cerr << "Error: Invalid number format\n";
return 1;
}
}
int n = 2;
if (argc > 2) {
n = atoi(argv[2]);
}
// Precision in bits (~3.3 bits per decimal digit)
PRECISION = static_cast<long>(DIGITS * 3.32193 + 16);
mpreal::set_default_prec(PRECISION); // Set default precision in bits
mpfr::mpreal::set_emin(-PRECISION - 2);
mpfr::mpreal::set_emax(PRECISION + 2);
mpreal x = 2.0;
if (argc > 1) {
x = mpreal(argv[1]);
}
auto start = chrono::high_resolution_clock::now();
mpreal nth_root = ::nth_root(x, n, DIGITS);
auto stop = chrono::high_resolution_clock::now();
assert(!isnan(nth_root));
chrono::duration<double> elapsed = stop - start;
cout << "x = " << x << endl;
cout << "nth_root(x, " << n << ") = " << nth_root.toString(DIGITS) << endl;
if (x > 0)
cout << "std::pow(x, 1/" << n << ") = " << setprecision(16) << std::pow(x.toDouble(), 1./n) << endl;
else if (n % 2 == 1)
cout << "std::pow(x, 1/" << n << ") = " << setprecision(16) << -std::pow(-x.toDouble(), 1./n) << endl;
printf("Computation Time: %f seconds\n", elapsed.count());
return 0;
}
% ./nth_root 2 12 100
Iterations: 215
x = 2
nth_root(x, 12) = 1.059463094359295264561825294946341700779204317494185628559208431458761646063255722383768376863945569
std::pow(x, 1/12) = 1.059463094359295
Computation Time: 0.000295 seconds
Computing \(ln x\)
We can also us the Newton-Raphson method to compute the natural logarithm. So to compute \(ln N\) we set \(f(x) = e^x - N\). The derivative is \(f’(x) = e^x\). Our iterative formula then becomes:
\[ x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)} = x_n - \frac{e^{x_n} - N}{e^{x_n}} = x_n - 1 + \frac{N}{e^{x_n}} \]
So to compute the next iteration of \(x\), we need to compute \(e^x\) to the desired number of digits. We can use the binary splitting version of exp(x) developed above or use the exp(x) in the mpreal library.
newton_ln.cpp
/**
* @file newton_ln.cpp
* @brief Program to compute ln x.
*
* This program can compute ln x using Newton–Raphson method.
*
* Syntax: newton_ln [x] [num_digits]
*
* get mpreal.h at https://github.com/advanpix/mpreal.git
*
* Compile with command:\n
* macOS:\n
* install homebrew https://brew.sh\n
* put mpreal.h in /opt/homebrew/include\n
* brew install gcc\n
* brew install gmp\n
* brew install mpfr\n
* g++ -std=c++20 -O3 -I/opt/homebrew/include -L/opt/homebrew/lib -o newton_ln newton_ln.cpp -lmpfr -lgmp
*
* Linux\n
* put mpreal.h in /urs/local/include\n
* sudo apt install libgmp-dev\n
* sudo apt install libmpfr-dev\n
* g++ -std=c++20 -O3 -o newton_ln newton_ln.cpp -lmpfr -lgmp
*
* @author Richard Lesh
* @date 2025-07-26
* @version 1.0
*/
#include <cassert>
#include <cmath>
#include <gmpxx.h> // GMP C++ Wrapper
#include <iostream>
#include <iomanip>
#include <mpreal.h> // MPFR C++ Wrapper
#include <ThreadPool.hpp>
#include <semaphore>
using namespace mpfr;
using namespace std;
long DIGITS = 50; // in decimal digits
long PRECISION; // in bits
long NUM_TERMS;
int MAX_THREAD_DEPTH = static_cast<int>(round(log(thread::hardware_concurrency()) / log(2.0)));
/**
* @brief Compute ln x.
*
* Function to compute \f( ln(x) \f) using Newton–Raphson method.
*
* @param x0 Argument to use when computing \f( ln(x) \f)
* @param digits Decimal digits of accuracy
* @return mpreal A multi-precision value for \f( ln(x) \f).
*/
// NOLINTNEXTLINE(misc-no-recursion)
mpreal ln_newton(const mpreal& x0, long num_digits) {
if (x0 <= 0) throw domain_error("ln(x) undefined for x <= 0");
// Reduce x0 to [1,2] range using: ln(x) = ln(m) + k * ln(2) where
// x = m * 2^k
int k = 0;
mpreal N = x0;
while (N > 2) {
N /= 2;
k++;
}
while (N < 1) {
N *= 2;
k--;
}
mpreal x = std::log(N.toDouble()); // reasonable initial guess
x.set_prec(PRECISION);
mpreal max_relative_err = pow(10, -num_digits);
mpreal ln2 = (k == 0) ? 0.0 : ln_newton(2.0, DIGITS);
int i = 0;
mpreal last;
do {
last = x;
x += (N / exp(x)) - 1;
++i;
} while (abs(x - last) / x > max_relative_err);
cout << "Iterations: " << i << endl;
return x + k * ln2;
}
/**
* @brief Main function.
*
* The entry point of the program.
*
* Syntax: newton_ln [x] [num_digits]
*
* @return int Returns 0 on successful execution.
*/
int main(int argc, char **argv) {
if (argc > 2) {
char* endptr;
DIGITS = static_cast<long>(strtod(argv[2], &endptr));
if (*endptr != '\0' || endptr == argv[2]) {
cerr << "Error: Invalid number format\n";
return 1;
}
}
// Precision in bits (~3.3 bits per decimal digit)
PRECISION = static_cast<long>(DIGITS * 3.32193 + 16);
mpreal::set_default_prec(PRECISION); // Set default precision in bits
mpfr::mpreal::set_emin(-PRECISION - 2);
mpfr::mpreal::set_emax(PRECISION + 2);
mpreal x = 2.0;
if (argc > 1) {
x = mpreal(argv[1]);
}
auto start = chrono::high_resolution_clock::now();
mpreal ln = ln_newton(x, DIGITS);
auto stop = chrono::high_resolution_clock::now();
assert(!isnan(ln));
chrono::duration<double> elapsed = stop - start;
cout << "x = " << x << endl;
cout << "ln(x) = " << ln.toString(DIGITS) << endl;
cout << "std::ln(x) = " << setprecision(16) << std::log(x.toDouble()) << endl;
printf("Computation Time: %f seconds\n", elapsed.count());
return 0;
}
Golden Ratio
The Golden Ratio, \(\phi\), is defined as the positive solution to the equation:
\[ \frac{a + b}{a} = \frac{a}{b} \]
Solving this gives:
\[ \phi = \frac{1 + \sqrt{5}}{2} \approx 1.6180339887\ldots \]
Using the \(\sqrt{N}\) function that we developed above, we can easily compute the Golden Ratio to any desired precision.
golden_ratio.cpp
/**
* @file golden_ratio.cpp
* @brief Program to compute the golden ratio.
*
* This program can compute the Golden Ratio using
* the Newton–Raphson method for square roots.
*
* Syntax: golden_ratio [num_digits]
*
* get mpreal.h at https://github.com/advanpix/mpreal.git
*
* Compile with command:\n
* macOS:\n
* install homebrew https://brew.sh\n
* put mpreal.h in /opt/homebrew/include\n
* brew install gcc\n
* brew install gmp\n
* brew install mpfr\n
* g++ -std=c++20 -O3 -I/opt/homebrew/include -L/opt/homebrew/lib -o golden_ratio golden_ratio.cpp -lmpfr -lgmp
*
* Linux\n
* put mpreal.h in /urs/local/include\n
* sudo apt install libgmp-dev\n
* sudo apt install libmpfr-dev\n
* g++ -std=c++20 -O3 -o golden_ratio golden_ratio.cpp -lmpfr -lgmp
*
* @author Richard Lesh
* @date 2025-08-07
* @version 1.0
*/
#include <cassert>
#include <cmath>
#include <gmpxx.h> // GMP C++ Wrapper
#include <iostream>
#include <iomanip>
#include <mpreal.h> // MPFR C++ Wrapper
using namespace mpfr;
using namespace std;
long DIGITS = 50; // in decimal digits
long PRECISION; // in bits
/**
* @brief Compute the Golden Ratio.
*
* Function to compute the Golden Ratio using Newton–Raphson for square roots.
*
* @param num_digits Accuracy in decimal digits
* @return mpreal A multi-precision value for \f( \phi \f).
*/
mpreal golden_ratio(long num_digits) {
mpreal ϕ = (1 + sqrt(5)) / 2;
return ϕ;
}
/**
* @brief Main function.
*
* The entry point of the program.
*
* Syntax: golden_ratio [num_digits]
*
* @return int Returns 0 on successful execution.
*/
int main(int argc, char **argv) {
if (argc > 1) {
char* endptr;
DIGITS = static_cast<long>(strtod(argv[1], &endptr));
if (*endptr != '\0' || endptr == argv[1]) {
cerr << "Error: Invalid number format\n";
return 1;
}
}
// Precision in bits (~3.3 bits per decimal digit)
PRECISION = static_cast<long>(DIGITS * 3.32193 + 16);
mpreal::set_default_prec(PRECISION); // Set default precision in bits
mpfr::mpreal::set_emin(-PRECISION - 2);
mpfr::mpreal::set_emax(PRECISION + 2);
auto start = chrono::high_resolution_clock::now();
mpreal ϕ = ::golden_ratio(DIGITS);
auto stop = chrono::high_resolution_clock::now();
assert(!isnan(ϕ));
chrono::duration<double> elapsed = stop - start;
cout << "ϕ = " << ϕ.toString(DIGITS) << endl;
printf("Computation Time: %f seconds\n", elapsed.count());
return 0;
}
See my HiDrive for files with digits of π and other math constants.