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" .

reference of functions in this page class reference
SLong functions SLong class
SDouble functions SDouble class
SComplex functions SComplex class
Other functions

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.