Multi-Precision Arithmetic by C++ with no use of assembler
SN library Copyright (C) 1999-2018 K.Tsuru
Reference of mathematical functions
Here I shall explain the fundamental usage of mathematical functions. Their
prototype declarations are given in "snmath.h" .
1. SLong functions
The argument(s) in the functions below can be "SInteger" class object, but a few cannot. They are commented as "(SLong only)".
combL
function "combL(n, k );" returns the combination number (binomial coefficients) nCk.
form SLong combL(unsigned long n, unsigned long k);
comment Give the numbers which satisfy the condition n>=k. It returns nCk. If n<k, an error occurs.
EulerNum
function "EulerNum(n);" returns nth Euler's number.
form SLong EulerNum(uint n );
gcdL
function "gcdL(a, b );" returns the greatest common divisor of two integers a and b .
form SLong gcdL(const SLong& a, const SLong& b);
Lpow
function "Lpow(x , n ); "returns x n.
form SLong Lpow(const SLong& x, unsigned long n);
comment
See a sample file "smppow2.cpp" for usage and Ipow() function below.
LRand
function "LRand(...);" generates a SLong random number.
form SLong LRand(unsigned int size, FType rdx = DRADIX);
comment Give the size to size and the radix (DRADIX or BRADIX) to rdx . It uses the standard functions "time()" and "rand()".
It always attaches the positive sign.
LSqrt(SLong only)
function "LSqrt(N);" returns the square root of a SLong integer N .
form SLong LSqrt(const SLong& N, int* sgn = NULL);
comment For a SLong integer N it returns the maximum integer q which satisfies the condition N-q2>=0. If you give the address to "sgn", you get the sign of N-q2. sgn>=0, when N is a perfect square integer sgn==0.
Fact
function "Fact(n )" returns the factorial of n i.e. n ! by the SLong integer.
form SLong Fact(unsigned long n);
MpowMod
function "MpowMod(M , d, n );" returns Md mod n .
form SLong MpowMod(const SLong& M, ulong d, const SLong& n);
Pow2
function "Pow2(n );" returns 2n.
form SLong Pow2(ulong n);
TanCoeffTable
function "TanCoeffTable(T, n);" makes a table of tangent coefficients in T up to nth.
form void TanCoeffTable(SNBlock <SLong>& T, uint n);
caution "Block" has been changed into "SNBlock" since ver.2.18.
comment See a sample program "smpbern.cpp" for usage.
2. SDouble functions
Acos
function "Acos(x );" returns an inverse trigonometric function arccos x .
form SDouble Acos(const SDouble& x );
comment For |x | <= 1.0 it returns the principle value of arccos x. When |x | > 1.0 an abnormal termination occurs in the Asin(x ) function used in it.
Mathematical constant values
Pi
function "Pi();" returns the value of pi.
form SDouble Pi();
comment AGM (Arithmetic-Geometric Mean) method is used. See the source file "sdfagmpi.cpp"
for algorithm. For other programs of pi see source files "sdc*pi.cpp"
or "sxc*pi.cpp".
Asin
function "Asin(x );" returns an inverse trigonometric function arcsin x .
form SDouble Asin(const SDouble& x );
comment For |x | <= 1.0 it returns the principle value of arcsin x. When |x | > 1.0 an abnormal termination occurs.
Algorithm of Asin
A. If x > 1/sqrt(2), using a formula
arcsin x = pi/2 - arcsin ( sqrt(1-x2) )
it can be reduced into |x | < 1/sqrt(2).
B. Pay attention a modified addition theorem
arcsin(x ) = arcsin(y ) + arcsin(delta ),
where x, y >=0 and
delta = x*sqrt(1-y2) - y*sqrt(1-x2).
Because this is an identity, the value of y can be chosen arbitrarily within the domain. Taking x = y gives delta = 0.
Here we think of the following idea.
1. In the double precision it evaluates x' = asin(x"). Where the x" in rhs is a rough estimated value of x in double precision if it is a multi-precision decimal number.
2. In the multi-precision it calculates y = Sin(x '). The convergence of series Sin(x') is fast. Here an equation
x' = arcsin(y)
stands up in the multi-precision.
3. It evaluates delta in the multi-precision.
4. Using a function AsinSeries(x) which provides the value of arcsin x by series,
Asin(x) = x' + AsinSeries(delta)
is calculated. Because the value of delta is very small less than 10-15, we get a fast convergence of series.
The algorithm mentioned above enables us to calculate arcsin x very rapidly. Most of processing time is spent on the evaluation Sin(x '). A similar method is used in the function Log(x ) and Atan(x ).
Atan
function "Atan(x );" returns the principle value of inverse trigonometric function
arctan x.
form SDouble Atan(const SDouble& x );
Atan2
function This is a SDouble version of atan2().
form SDouble Atan2(const SDouble& y, const SDouble& x);
comment It returns arctan(y /x ) whose value is between -pi and pi.
Acosh
function "Acosh(x ):" returns an inverse hyperbolic function arcosh x .
form SDouble Acosh(const SDouble& x);
Asinh
function "Asinh(x );" returns an inverse hyperbolic function arsinh x .
form SDouble Asinh(const SDouble& x);
Atanh
function "Atanh(x );" returns an inverse hyperbolic function artanh x .
form SDouble Atanh(const SDouble& x);
Cbrt
function "Cbrt(x );" returns the cubic root of x .
form SDouble Cbrt(const SDouble& x);
comment It returns x*RecCbrt(x*x).
Cos
function "Cos(x );" returns a trigonometric function cos x .
form SDouble Cos(const SDouble& x);
<new>CosBS since ver. 2.18
function "CosBS(x );" returns a trigonometric function cos x .
form SDouble CosBS(const SDouble& x);
see CosSinBS, ExpBS, SinBS, TanBS
Dpow
function "Dpow(x , n );" returns xn, where n is an integer and may be negative.
form SDouble Dpow(const SDouble& x, long n);
Dpow10
function"Dpow(n );" returns 10n, where n is an integer and may be negative.
form SDouble Dpow10(long n);
DpowD
function "DpowD(x , p );" returns x p (x >0), where p is a multi-precision real number.
form SDouble DpowD(const SDouble& x, const SDouble& p);
caution When x <=0 an error occurs. The multi-precision complex is not prepared.
DRand
function "DRand(...);" generates a SDouble random number.
form SDouble DRand(uint size, FType rdx = DRADIX, int exp = 0);
comment Give the size to size, the radix (DRADIX or BRADIX) to rdx and the exponent to exp. It uses the standard functions "time()" and "rand()".
It always attaches the positive sign.
E(constant)
function It provides the base of natural logarithm e in a decimal system.
form SDouble E(const SDouble* e = NULL, SDouble (*pfCalcFunc)() = SNE);
comment
If you have e which has enough significant figures, give it to e . It can be given from the file using SetConstByFile()below. If you have a function which is faster than that of this library,
give the pointer of that function to pfCalcFunc . In default it calls SNE(). Other constant functions Pi and Log10(contant) have the same form. Usually you can use these constant functions with no argument, such as
E(), Pi() and Log10().
Exp
function "Exp(x );" returns an exponential function exp(x).
form SDouble Exp(const SDouble& x );
comment It returns "ExpDS(x)" for rather short decimal and "ExpL(x)" for long one in below.
<new>ExpBS since ver. 2.18
form SDouble ExpBS(const SDouble& x);
comment It evaluates exp(x) using the binary splitting method.
ExpDS
function "ExpDS(x);" returns exp(x).
form SDouble ExpDS(const SDouble& x);
comment It uses a divide and series method. See comment in the source file "sdfexpx.cpp"
for details.
ExpL
function "ExpL(x);" returns exp(x).
form SDouble ExpL(const SDouble& x);
Algorithm of ExpL
Pay attention to an identity
exp(x) = r exp{ x - log(r) }.
1. It evaluates an approximate value r = ExpDS(x') in f figures.
2. Next a = x - log(r) is evaluated in full precision (F figures). In the evaluation log(x) for rather long decimal x the LogAGM() function (see below) is used.
a has an order of R-f where R is radix..
3. exp(a) is evaluated by series b = ExpSeries(a) in which the number of necessary terms is an order of F/f.
4. It returns br.
Expower
form SDouble Expower(long n);
comment "Expower(n );" returns en , where n is an integer and may be negative.
Fmod
function This is a SDouble version of fmod().
form SDouble Fmod(const SDouble& x, const SDouble& y, SDouble& a);
comment
It calculates a and f which satisfy x = ay + f (a is an integer and |f |<|y |). It returns f and sets a in third argument.
usage
SDouble x, y, f, a;
...............
f = Fmod(x, y, a); //not "&a".
Log
function "Log(x );" returns a natural logarithm log x .
form SDouble Log(const SDouble& x);
comment It returns "LogE(x)" for rather short decimal and "LogAGM(x)"
for long one in below.
LogE
function "LogE(x );" returns a natural logarithm log x .
form SDouble LogE(const SDouble& x);
LogAGM
function "LogAGM(x );" returns a natural logarithm log x for small x (< 1.0).
form SDouble LogAGM(const SDouble& x);
commentAGM (Arithmetic-Geometric Mean) method is used. It needs the constant Pi().
See the source file "sdflogag.cpp" for algorithm. It does not
check the range of value x. When x has a large value, an overflow error will occur.
Algorithm of LogE
It is well known that the convergence of logarithmic series is very slow.
It takes much time in a normal method. Then I introduce the following idea.
Paying attention to an identity
x = exp(r ){1-(1-x exp(-r))} ,
we obtain
log x = r + log(1-d)
where d= 1- x exp(-r).
The value of r can be arbitrarily chosen because of identity. The solution of an equation
d=0 is r= log x, then evaluating it in double precision the d has a very small value, order of 10-15.
Further solving an equation (1+u)/(1-u) = 1-d with respect to u , we obtain
u = -d/(2.0-d);
and
log(1-d) = log{(1+u)/(1-u)}= 2(u + u 3/3 + u5/5 +...).
The square of u has an order of 10-30 for any x (or r), then in order to obtain the value in 400 figures it is necessary to
sum up only about 15 terms.
In the practical program an intermediate step is inserted in which the
approximate value of r is evaluated in rather large figures. The processing speed depends on
that of exponential function Exp(x). When the value of x is too large or too small to treat by double precision, the processing
time of Log(10) is added if it has not been evaluated yet. In other case
it does not depend on the value of x .
Log10(constant)
function "Log10();" returns the natural logarithm of ten log(10)
form SDouble Log10(const SDouble* ln10=NULL, SDouble (*pfCalcFunc)()=SNLog10);
Log10
function "Log10(x );" returns the common logarithm log10 x .
form SDouble Log10(const SDouble& x );
Modf
function This is a SDouble version of modf().
form SDouble Modf(const SDouble& x, SDouble& ipart);
comment It sprits x into an integral part ipart and a decimal part d. It returns d and sets ipart in second argument.
usage
SDouble d, x, ipart;
...............
d = Modf(x, ipart); //not "&ipart".
Qrrt
function "Qrrt(x );" returns the quartic root of x (x>=0).
form SDouble Qrrt(const SDouble& x);
comment It returns x*{ RecQrrt(x) }3.
Sin
function "Sin(x );" returns a trigonometric function sin x .
form SDouble Asin(const SDouble& x);
Tan
function "Tan(x );" returns a trigonometric function tan x .
form SDouble Atan(const SDouble& x);
Cosh
function "Cosh(x );" returns a hyperbolic function cosh x .
form SDouble Cosh(const SDouble& x);
Sinh
function "Sinh(x );" returns a hyperbolic function sinh x .
form SDouble Sinh(const SDouble& x);
Tanh
function "Tanh(x );" returns a hyperbolic function tanh x .
form SDouble Tanh(const SDouble& x);
SetConstByFile
function It sets a constant on memory by reading from a text file.
form
void SetConstByFile(const char* fname, SDouble (*pfunc)(const SDouble* c, SDouble (*pf)()), SDouble (*pfCalcFunc)(), uint ef);
comment
1. Give a file name to fname.
2. Set pfunc = Log10, E or Pi. It does not check that the given value is correct or not.
3. Set pfCalcFunc a function pointer for the calculation. If a sufficient length (accuracy)
can been read from the file, NULL is ok. If pfCalcFunc == NULL and the length is not enough, a syntax error occurs in EntryConst().
4. Set ef your significant figures in DRADIX. If ef = 0, it sets current significant figures given by SetEffFig() .
example
SetConstByFile("pi.snc", Pi, NULL, 10000u);
10000*4 = 40000 digits in decimal are read from "pi.snc" file and set in Pi().
RecCbrt
function "RecCbrt(x );" returns x -1/3.
form SDouble RecCbrt(const SDouble& x);
RecQrrt
function "RecQrrt(x );" returns x -1/4(x>0).
form SDouble RecQrrt(const SDouble& x);
RecSqrt
function "RecSqrt(x ); " returns 1/sqrt(x ) (x>0).
form SDouble RecSqrt(const SDouble& x );
commentIt provides 1/sqrt(x) using Newton's method i.e. by solving an equation
x - 1/y2 = 0.
Algorithm of RecSqrt
Taking w=x *100pthe value of p is chosen to satisfy a condition 0.01<w <1. The initial value is taken as
y0= 1/sqrt(w) (calculated by double).
It repeats
yn+1= yn + dy where dy= 0.5yn(1-w yn2).
Here y has a value y = 1/sqrt(x*100p)= 1/(sqrt(x)*10p), then considering 1/sqrt(x ) = y*10p, it multiplies y by 10p at last.
If the "doubling significant figures" method is used, it becomes a few times faster. It is very efficient as
yn = 0.abcd ..... efgh ijkl
dy = 0.0000 ..... 0000 EFGH IJKL ....,
the non-zero part of dy overlaps with that of yn by a few figures.
Sqrt
function "Sqrt(x );" returns the square root of SDouble x (x>=0)
form SDouble Sqrt(const SDouble& x);
comment It returns x*RecSqrt(x ).
Pi(constant)
function It returns the value of pi in the current significant figures.
form SDouble Pi(const SDouble* pi=NULL, SDouble (*pfCalcFunc)() = SNPi);
3. SComplex functions
The SComplex class provides a multi-precision complex number arithmetic.
I referred to gcc's "complext.h".
Arg
function "Arg(z ); " returns the argument of z in xy plane by a SDouble value.
form SDouble Arg(const SComplex& z );
Cabs
function"Cabs(z ); " returns the absolute value of z i.e. |z| by a SDouble value.
form SComplex Cabs(const SComplex& z );
Cexp
function"Cexp(z ); " returns exp(z ) .
form SComplex Cexp(const SComplex& z );
Clog
function "Clog(z ); " returns log(z ) .
form SComplex Clog(const SComplex& z );
Cpow
function "Cpow(x, y ); " returns x y.
form
SComplex Cpow(const SComplex& z, int n);
SComplex Cpow(const SComplex& x, const SDouble& y);
SComplex Cpow(const SDouble& x, const SComplex& y);
SComplex Cpow(const SComplex& x, const SComplex& y);
Conj
function "Conj(z); " or "z.Conj();" returns the complex conjugate of z.
form SComplex Conj(const SComplex& z);
SComplex SComplex::Conj()
const;
comment Let z = x + i y it returns z *= x - i y.
Norm
function "Norm(z); " or "z.Norm();" returns the square of the magnitude of z. by a SDouble value.
form SDouble Norm(const SComplex& z);
SDouble SComplex::Norm()
const;
comment Let z = x + i y, it evaluates x2 + y2.
Polar
function "Polar(r, a ); " returns r{ cos(a) + i sin(a) }.
form SComplex Polar(const SDouble& r,const SDouble& a);
Real
function "Real(z ); " or "z.Real();" returns the real part of z by a SDouble value.
form SDouble Real(const SComplex& z);
inline SDouble SComplex::Real()
const;
Imag
function "Imagl(z); " or "z.Imag();"returns the imaginary part of z by a SDouble value.
form SDouble Imag(const SComplex& z);
inline SDouble SComplex::Imag()
const;
ImagUnit(constant)
function "ImagUnit" returns the imaginary unit i.e. sqrt(-1)
= SComplex(0, 1).
form extern const SComplex ImagUnit;
Csqrt
function "Csqrt(z ); " returns the square root of z.
form SComplex Csqrt(const SComplex& z);
Ccos
function "Ccos(z ); " returns cos(z ) .
form SComplex Ccos(const SComplex& z);
Csin
function "Csin(z ); " returns sin(z ) .
form SComplex Csin(const SComplex& z);
Ctan
function "Ctan(z ); " returns tan(z ) .
form SComplex Ctan(const SComplex& z);
Ccosh
function "Ccosh(z ); " returns cosh(z ) .
form SComplex Ccosh(const SComplex& z);
Csinh
function "Csinh(z ); " returns sinh(z ) .
form SComplex Csinh(const SComplex& z);
Ctan
function "Ctanh(z ); " returns tanh(z ) .
form SComplex Ctanh(const SComplex& z);
4. Other functions
BernNum
function "BernNum(n );" returns n-th Bernouilli's number by a SFraction value.
form SFraction BernNum(uint n );
Ipow
function "Ipow(x , n ); "returns x n.
form SInteger Ipow(const SInteger& x, unsigned long n );
comment
This is a SInteger version of Lpow(). See a sample file "smppow2.cpp" for usage.