Numerical Methods at work

Disclaimer:
Permission to use, copy, and distribute this software, and It’s documentation for any non-commercial purpose is hereby granted without fee, provided: THE SOFTWARE IS PROVIDED "AS-IS" AND WITHOUT WARRANTY OF ANY KIND, EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. IN NO EVENT SHALL Henrik Vestermark, BE LIABLE FOR ANY SPECIAL, INCIDENTAL, INDIRECT OR CONSEQUENTIAL DAMAGES OF ANY KIND, OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER ADVISED OF THE POSSIBILITY OF DAMAGE, AND ON ANY THEORY OF LIABILITY, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.

Arbitrary precision package. (Minor Release July 2023)

Arbitrary precision for integers, floating points, complex numbers, etc. Nearly everything is here! A collection of five C++ header files. The news is that the internal format has changed from storing the number as a string to storing the number as binary digits. This has resulted in a huge performance gain. Furthermore, in this maintenance release (September 2022) I have added new algorithms for sqrt(), inverse, π, and e and also improved all the trigonometric and hyperbolic functions, plus conversion to and from a string, representation has seen significant performance gains and lastly, many minor bugs has been fixed. In the previous January 2023 release I have added all the specials functions and constants (see below). For this release in July 2023 a random precision class has been added.

  • Arbitrary integer precision class
  • Arbitrary Floating-point precision class
  • A portable Complex numeric template<class T>
  • A portable interval arithmetic template<class T>
  • A portable fraction arithmetic template<class T>
Functions & Operators for the arbitrary precision class.
  • All standard C++ operators: +,-,*,/,%,++(both pre and postfix),--(both pre and postfix),[], ≪, ≫
  • Logical operators: &, |, ^, &&, ||, !
  • Assignment operators: =,+=, -=, *=, /=, %=, ≪=, ≫=, &=, |=, ^=
  • Comparison operators: ==, !=, ≤, ≥, <, >
  • Logarithm & Exponential functions: exp(), log(), log10()
  • Trigonometric functions: sin(), cos(), tan(), atan(), asin(), acos(), atan2()
  • Miscellaneous: pow() & sqrt(), n-rooth of a number, gcd() Greatest coming divisor, lcm() least common multiple
  • Hyperbolic functions: sinh(), cosh(), tanh(), asinh(), acosh(), atanh()
  • Special functions: tgamma(), beta(), erf(), erfc(), zeta(), Bernoulli numbers, Bernoulli polynomials, factorials, falling factorials & binomials
  • Universal Arbitrary precisions constants π,√2, √3 ,1/√2, 1/√3, LN2, LN10, e, Euler Mascheroni γ, Catalan and Apéry constant ζ(3)
  • Supporting functions for formatting floating-point output: toString(), toFixed(), toPrecision() π toExponential() (borrow from JavaScript methods)
  • A random precision integer class. You can now generate arbitrary precision random numbers using well-known random classes like Mersenne Twister, Xoshiro 256** and 512**, plus the cryptographic graded ChaCha20 and others
  • Several other methods (see manual)

Image of infinite digits

Furthermore, for each floating precision number, the working rounding mode for arithmetic operations can be controlled.
Four rounding modes are supported:

The 4 rounding modes make it easy to implement interval arithmetic, which means you can now get a precise bound of the error for every floating point calculation.
Technically the number of digits for a number that can be handled is around 4 Billions digits, however, most likely you will run into system limitations before that. However, I have been working with numbers that exceed 100-500 million digits without any issues!
Click here to download the user manual Download
Also, don't forget to check out our document on the math behind arbitrary precision. Click here for Download


Why use this package instead of Gnu's GMP?

  • It has less restrictive permission rules.
  • It supports all relevant trigonometric, logarithms and exponential  functions like exponential, logarithm, trigonometric and hyperbolic functions, etc. which GMP does not
  • It's born as a C++ class  and not a C library with a C++ wrapper.
  • You also have rounding controls and can freely be mixed among variables with different precisions
  • π, √2, √3 ,1/√2, 1/√3, Ln2, Ln10 is available in arbitrary precision.
  • It is multithreaded for improved performance
  • Easier to use

Why use Gnu's GMP?

  • Because it's GNU!
  • Faster and more choices on basic functions and algorithms
  • Gnu's GMP can be located at: www.gnu.org/software/gmp

Please note that I did not develop this package to compete with Gnu's GMP but rather because I was missing features not found in GMP.
However, since I get a lot of questions why? I have tried to answer it above. Have fun.

Download source file for the newest version

Arbitrary Precision Math C++ Package.  Include C++ header file & pdf documentation.

The current version is: (updated july 2023)

In this maintenance release of the new "Binary" version, it has been properly tested in a 64-bit environment with Visual Studio 2022 C++ compiler and with Code::Blocks 20.03 GNU Compiler. All GNU compiler warnings and errors have been fixed for better portability. Many critical functions have been optimized for better performance.
In the previous release, the sqrt() is 3 times faster using a new algorithm. the n-th root of a number is now available as a special function nroot(). Calculation of π is now many times faster using the Chudnovsky binary splitting algorithm and the transcendental constants ln(2), ln(10), exp(1) is now 20-100 times faster than in the previous version using a spigot algorithm or binary splitting method. New algorithms have been developed to speed up the conversion from string to the internal binary format and vice versa from binary format to string. Lastly, all the trigonometric and hyperbolic functions have added coefficient scaling.(see the papers related to these new methods on the left side menu).

A special thanks to Robert McInnes who help improve and test the performance of the π implementation and was instrumental to test the accuracy of π to more than 600 million digits.


Examples of using int_precision & float_precision:

{
int_precision a(11111111), d;
float_precision f1(10,64); // Initialized to zero and set precision to 64 digits
float_precision f2(2.0), f3(3), f4("0.12345678901234567890" );

d = a * a * a; // a^3
cout << "11111111^3=" << d << endl;
d /= a;
cout << "11111111^2=" << d << endl;

f1 = log( f1 );
f1 *= f2 + f3 / f4;
f1 = _float_table( _PI, 48 ); // Get PI with 48 decimals
cout << "PI with 48 digits:" << f1 << endl;
f1 = float_precision( 2 );
f1 = sqrt( f1 );
cout << "Sqrt(2) with 48 digits:" << f1 << endl;
}

Miscellaneous functions

For integer precision the following extra functions are available:

  • int_precision ipow( const int_precision&, const int_precision& ) // a^b
  • int_precision ipow_modular( const int_precision&, const int_precision&, const int_precision& ) // a^b%c
  • bool iprime( const int_precision& ) // Fast function to check if a int_precision number is a prime.
  • gcd(const int_precision&, const int_precision&) // Greatest Common Divisor
  • lcm(const int_precision&, const int_precision&) // Least Common multiple

For floating point precision the following extra functions are available (equivalent to the same function in the standard C library):

  • float_precision modf( float_precision, float_precision * )
  • float_precision fmod( float_precision, float_precision )
  • float_precision floor( float_precision )
  • float_precision ceil( float_precision )
  • float_precision frexp( float_precision, int * )
  • float_precision ldexp( float_precision, int )

Example using complex_precision and float_precision

// This example shows how to evaluate a polynomial with complex coefficients an at a complex point z. The evaluation is done using
// the precision of z and using the Horner algorithm

#include "iprecision.h"
#include "fprecision.h"
#include "complexprecision.h"

complex_precision<float_precision> horner( const register int n, const complex_precision<float_precision> a[],const complex_precision<float_precision> z )
{ complex_precision<float_precision> fval(0,z); // fval will inherit the same precision as z
fval = a[ 0 ];
for( int i = 1; i <= n; i++ )
fval = fval * z + a[ i ];

return fval;
}

Example The Brent-Salamin algorithm for calculating π to arbitrary precision

// Begin Brent-Salamin
std::string pi_brent_salamin(unsigned int digits)
{
const unsigned int extra = 5;
const int limit = -(int)(digits+2);
const float_precision c0(0), c1(1), c2(2), c05(0.5);
float_precision a(1, digits+extra), b(2,digits+extra), sum(0.5,digits+extra);
float_precision ak(0, digits+extra), bk(0, digits+extra), ck(1, digits+extra);
float_precision ab(0, digits+extra), asq(0, digits+extra);
float_precision pow2(1,digits), pi(3, digits+extra);
std::string ss;
b = c1 / sqrt(b);
for ( ; ck != c0 && ck.exponent()>limit; )
{

ak = c05*(a + b);
ab = a * b;
bk = sqrt(ab);
asq = ak * ak;
ck = asq - ab;
pow2 *= c2;
sum -= pow2*ck;
a = ak; b = bk;
}

pi = c2 * asq / sum;
pi.precision(digits);
ss = pi.toString();
return ss;
}

// End Brent-Salamin

There is also a small web tool where you can see how different multiplication algorithm works for arbitrary precision using the author Javascript BigFloat library. It supports a number of different methods, see below:

  • Karatsuba (2-way splitting).
  • Toom-cook3 (3-way splitting).
  • Schönhage-Strassen with Linear convolution.
  • Fast Fourier Multiplication.
  • School Book multiplication.
Click Here for the multiplication calculator