- /* snmath.h by K.Tsuru */
- /*********************
- SN library
- Mathematical functions since version 2.30
- Separate from "snfunc.h" for convenience.
- Almost trigonometric and exponential functions are supplied by the binary splitting,
- the E.A.Karatuba's method and the Brent's trick. For the trigonometric functions the "SComplex" or
- "SLComplex" arithmertic is applied to exp(i*x).
- **********************/
- #ifndef SN_MATH_H
- #define SN_MATH_H
- #ifndef SN_FUNCTIONS_H
- #include "snfunc.h"
- #endif // SN_FUNCTIONS_H
- ////////////// SLong class /////////////////
- // greatest common divisor
- SLong gcdL(const SLong& a, const SLong& b); // 2000
- SLong Lpow(const SLong& x, ulong n); // n-th power of x(x^n) 2002
- //10^n
- SLong Lpow10(long n); // "friend" is deleted since version 2.192
- // make a SLong random number in radix rdx
- SLong LRand(uint size, FType rdx = DRADIX); // 2003
- // square root of SLong integer N, q = [sqrt(N)]
- // sgn has sign of N - q*q, then sgn >= 0
- // If N is a square of an integer, sgn = 0.
- SLong LSqrt(const SLong& N, int* sgn = NULL); // 2004
- // return M^d mod n
- SLong MpowMod(const SLong& M, ulong d, const SLong& n);// 2005
- SLong Pow2(ulong p); // 2^p 2006
- // If m and n have a common divisor in type 2^b*R^r, they are divided by it.
- void ReduceByPow2Rdx(SLong& m, SLong& n, SLong* divisor = NULL);// 2007
- // least common multiple
- SLong lcmL(const SLong& a, const SLong& b); // 2008 since ver 2.8
- inline SLong Fact(int n) { return FactPF(n); } // factorial n! since version 2.21
- inline SLong binomial(ulong n, ulong r) { return combPF(n, r); } // same as below
- inline SLong comb(ulong n, ulong k) {return combPF(n,k); } // combination number(binomial) nCk
- inline SLong perm(const ulong n, const ulong r) { return permL(n, r); } // permutation number nPr
- SLong EulerNum(uint n); // n-th Euler number 4101
- SLong TanCoeff(uint n); // n-th tangent coefficient 4104
- void TanCoeffTable(SNBlock <SLong>& T, uint N); //table of tangent coefficient 4105
- /////////////////// SFraction class /////////////////
- SFraction BernNum(uint n); // Bernouilli's number 7100
- void BernNumTable(SNBlock <SFraction>& Bn, uint N); // table of Bernouilli's number 7101
- // Digit figures of factorial n! by Stirling's formula
- unsigned long Stirling(long n); // ver. 2.17.1 change into "unsigned long" from "long".
- /////////////////// SDouble class /////////////
- // Convert SDouble to an approximate in double.
- // If set err = 0, return DBL_MAX for overflow or zero for underflow error in double.
- // Error ID(OVERFLOW_ERR or UNDERFLOW_ERR) can be obtained by SNManeger::SNError().
- double doubleD(const SDouble& sd, int err = 1);// convert to double 3001
- SDouble Dpow(const SDouble& x, long n); //p-th power of x(x^p) 3002
- SDouble DpowD(const SDouble& x, const SDouble& p); //3003
- SDouble Dpow10(long n); // 10^n 3004
- // make a SDouble random number in radix rdx
- SDouble DRand(uint size, FType rdx = DRADIX, int exp = 0);// 3005
- inline SDouble Drand(uint size, FType rdx = DRADIX, int exp = 0) { // small capital version
- return DRand(size, rdx, exp);
- }
- /*-----------------
- divide x into integer and decimal part
- SDouble version of modf()
- x = ipart + dpart
- For example
- x = -1.5
- ipart = -1.0
- It returns -1.5 - (-1) = -0.5.
- --------------------*/
- SDouble Modf(const SDouble& x, SDouble& ipart); // 3008
- // SDouble version of fmod()
- SDouble Fmod(const SDouble& x, const SDouble& y, SDouble& ipart); // 3009
- // square root
- SDouble RecSqrt(const SDouble& x); // 1/sqrt(x) 3006
- SDouble Sqrt(const SDouble& x); // sqrt(x) 3007
- // cubic root
- SDouble Cbrt(const SDouble& x); // x^(1/3) 3010
- // reciporocal cubic root
- SDouble RecCbrt(const SDouble& x); // x^(-1/3) 3011
- // quartic root
- SDouble Qrrt(const SDouble& a); // x^(1/4) 3012
- // reciporocal quartic root
- SDouble RecQrrt(const SDouble& a); // x^(-1/4) 3013
- // exponential and hyperbolic functions
- SDouble Expower(long n); // n-th power of e 3302
- /**************************************************************************
- trigonometric functions calculated by DRAIDX
- **************************************************************************/
- SDouble Acos(const SDouble& x); // 3101
- SDouble Asin(const SDouble& x); // 3102
- inline SDouble Atan(const SDouble& x) { return AtanT(x); } // since version 2.21
- inline SDouble Cos(const SDouble& x) { return CosBS(x); } // since version 2.21
- inline SDouble Sin(const SDouble& x) { return SinBS(x); }
- inline SDouble Tan(const SDouble& x) { return TanBS(x); }
- // SDouble version of C's "double atan2(double y, double x)" arctan(y/x) 3110
- SDouble Atan2(const SDouble& y, const SDouble& x);
- /**************************************************************************
- It provides the exponential function exp(x).
- From the condition "exp(x) < DRADIX^DRADIX_EXP_MAX" the range of 'x' must be
- within
- |x| < DFIGURES*log(10)*DRADIX_EXP_MAX = 301759.17...
- to inline since version 2.21
- ***************************************************************************/
- inline SDouble Exp(const SDouble& x) {
- return ExpBS(x);
- }
- inline SDouble Cosh(const SDouble& x) {
- return CoshBS(x);
- }
- inline SDouble Sinh(const SDouble& x) {
- return SinhBS(x);
- }
- inline SDouble Tanh(const SDouble& x){
- return TanhBS(x);
- }
- // inverse hyperbolic functions
- SDouble Acosh(const SDouble& x); // 3309
- SDouble Asinh(const SDouble& x); // 3310
- SDouble Atanh(const SDouble& x); // 3311
- // logalithm functions
- inline SDouble Log(const SDouble& x) { return LogE(x); } // log(x)
- inline SDouble Log10(const SDouble& x){ return Log(x)*RecLog10(); } // log_10(x)
- /***********************************************************************
- Prototype declaration of constant functions has a formula
- SDouble func(const SDouble* userV = NULL,SDouble (*pfCalcFunc)() = fname);
- userV : a value given by user
- pfCalcFunc : default function, NULL is ok. A funcion which is faster than that
- of SN library can be used.
- See also "SetConstByFile()" in below.
- ************************************************************************/
- SDouble E(const SDouble* e = NULL,
- SDouble (*pfCalcFunc)(void) = SNE); // base of natural logarithm 3521
- SDouble Log10(const SDouble* ln10 = NULL,
- SDouble (*pfCalcFunc)(void) = SNLog10); // log(10.0) 3522
- SDouble Pi(const SDouble* pi = NULL,
- SDouble (*pfCalcFunc)(void) = SNPi); // pi 3523
- SDouble M2Pi(void); // 2/pi 3506
- SDouble RecE(void); // 1/e 3507
- SDouble MPi2(void); // pi/2 3511
- SDouble MPi4(void); // pi/4 3550
- SDouble RecLog10(void); // 1/log(10.0) 3551
- /////////////////// SComplex class /////////////////
- #ifndef SCOMPLEX_H
- #include "scomplex.h"
- #endif // SCOMPLEX_H
- inline SComplex Conj(const SComplex& z) { return z.Conj(); } // to inline since version 2.21
- inline SDouble Cabs(const SComplex& z) { return Sqrt(z.Norm()); } // |z| most simple
- inline SDouble Arg(const SComplex& z) { return Atan2(z.Imag(), z.Real()); } // argument in xy plane
- inline SDouble Real(const SComplex& z){ return z.Real();} // real part
- inline SDouble Imag(const SComplex& z){ return z.Imag();} // imaginary part
- inline SDouble Re(const SComplex& z) { return z.Real();} // Re(z) real part
- inline SDouble Im(const SComplex& z) { return z.Imag();} // Im(z) imaginary part
- inline SDouble Norm(const SComplex& z){ return z.Norm(); } // square of the magnitude
- /*********************************
- SComplex mathematical functions
- Acooding to GCC's "complex.h" the initial letter
- of funtion name is 'C';
- **********************************/
- SComplex Csqrt(const SComplex& z); // sqrt(z) 9101
- SComplex Cpolar(const SDouble& r, const SDouble& angle); // r*{cos(angle)+i*sin(angle)}; 9102
- SComplex Cpow(const SComplex& z, int n); // z^n 9103
- SComplex Cpow(const SComplex& x, const SDouble& y);// x^y 9104
- SComplex Cpow(const SDouble& x, const SComplex& y);// x^y 9105
- SComplex Cpow(const SComplex& x, const SComplex& y);// x^y 9106
- SComplex Cexp(const SComplex& z); // exp z 9107
- SComplex Clog(const SComplex& z); // log z 9108
- SComplex Ccosh(const SComplex& z); // cosh z 9109
- SComplex Csinh(const SComplex& z); // sinh z 9110
- SComplex Ctanh(const SComplex& z); // tanh z 9111
- SComplex Ccos(const SComplex& z); // cos z 9112
- SComplex Csin(const SComplex& z); // sin z 9113
- SComplex Ctan(const SComplex& z); // tan z 9114
- //The following are added since ver. 2.18.
- SComplex Cacos(const SComplex& z); // arccos z 9115
- SComplex Casin(const SComplex& z); // arcsin z 9116
- SComplex Catan(const SComplex& z); // arctan z 9117
- SComplex Cacosh(const SComplex& z); // arcosh z 9118
- SComplex Casinh(const SComplex& z); // arsinh z 9119
- SComplex Catanh(const SComplex& z); // artanh z 9120
- // See a sample program "smpcmpx3.cpp".
- // complex <double> complexC(const SComplex& z); //9151 ver. 2.17
- // conversion SComplex to complex <double>
- inline complex <double> complexC(const SComplex& z){
- return complex <double>(doubleD(z.Real()), doubleD(z.Imag()));
- }
- #endif // SN_MATH_H