1. /* snmath.h by K.Tsuru */
  2. /*********************
  3. SN library
  4. Mathematical functions since version 2.30
  5. Separate from "snfunc.h" for convenience.
  6. Almost trigonometric and exponential functions are supplied by the binary splitting,
  7. the E.A.Karatuba's method and the Brent's trick. For the trigonometric functions the "SComplex" or
  8. "SLComplex" arithmertic is applied to exp(i*x).
  9. **********************/
  10. #ifndef SN_MATH_H
  11. #define SN_MATH_H
  12. #ifndef SN_FUNCTIONS_H
  13. #include "snfunc.h"
  14. #endif // SN_FUNCTIONS_H
  15. ////////////// SLong class /////////////////
  16. // greatest common divisor
  17. SLong gcdL(const SLong& a, const SLong& b); // 2000
  18. SLong Lpow(const SLong& x, ulong n); // n-th power of x(x^n) 2002
  19. //10^n
  20. SLong Lpow10(long n); // "friend" is deleted since version 2.192
  21. // make a SLong random number in radix rdx
  22. SLong LRand(uint size, FType rdx = DRADIX); // 2003
  23. // square root of SLong integer N, q = [sqrt(N)]
  24. // sgn has sign of N - q*q, then sgn >= 0
  25. // If N is a square of an integer, sgn = 0.
  26. SLong LSqrt(const SLong& N, int* sgn = NULL); // 2004
  27. // return M^d mod n
  28. SLong MpowMod(const SLong& M, ulong d, const SLong& n);// 2005
  29. SLong Pow2(ulong p); // 2^p 2006
  30. // If m and n have a common divisor in type 2^b*R^r, they are divided by it.
  31. void ReduceByPow2Rdx(SLong& m, SLong& n, SLong* divisor = NULL);// 2007
  32. // least common multiple
  33. SLong lcmL(const SLong& a, const SLong& b); // 2008 since ver 2.8
  34. inline SLong Fact(int n) { return FactPF(n); } // factorial n! since version 2.21
  35. inline SLong binomial(ulong n, ulong r) { return combPF(n, r); } // same as below
  36. inline SLong comb(ulong n, ulong k) {return combPF(n,k); } // combination number(binomial) nCk
  37. inline SLong perm(const ulong n, const ulong r) { return permL(n, r); } // permutation number nPr
  38. SLong EulerNum(uint n); // n-th Euler number 4101
  39. SLong TanCoeff(uint n); // n-th tangent coefficient 4104
  40. void TanCoeffTable(SNBlock <SLong>& T, uint N); //table of tangent coefficient 4105
  41. /////////////////// SFraction class /////////////////
  42. SFraction BernNum(uint n); // Bernouilli's number 7100
  43. void BernNumTable(SNBlock <SFraction>& Bn, uint N); // table of Bernouilli's number 7101
  44. // Digit figures of factorial n! by Stirling's formula
  45. unsigned long Stirling(long n); // ver. 2.17.1 change into "unsigned long" from "long".
  46. /////////////////// SDouble class /////////////
  47. // Convert SDouble to an approximate in double.
  48. // If set err = 0, return DBL_MAX for overflow or zero for underflow error in double.
  49. // Error ID(OVERFLOW_ERR or UNDERFLOW_ERR) can be obtained by SNManeger::SNError().
  50. double doubleD(const SDouble& sd, int err = 1);// convert to double 3001
  51. SDouble Dpow(const SDouble& x, long n); //p-th power of x(x^p) 3002
  52. SDouble DpowD(const SDouble& x, const SDouble& p); //3003
  53. SDouble Dpow10(long n); // 10^n 3004
  54. // make a SDouble random number in radix rdx
  55. SDouble DRand(uint size, FType rdx = DRADIX, int exp = 0);// 3005
  56. inline SDouble Drand(uint size, FType rdx = DRADIX, int exp = 0) { // small capital version
  57. return DRand(size, rdx, exp);
  58. }
  59. /*-----------------
  60. divide x into integer and decimal part
  61. SDouble version of modf()
  62. x = ipart + dpart
  63. For example
  64. x = -1.5
  65. ipart = -1.0
  66. It returns -1.5 - (-1) = -0.5.
  67. --------------------*/
  68. SDouble Modf(const SDouble& x, SDouble& ipart); // 3008
  69. // SDouble version of fmod()
  70. SDouble Fmod(const SDouble& x, const SDouble& y, SDouble& ipart); // 3009
  71. // square root
  72. SDouble RecSqrt(const SDouble& x); // 1/sqrt(x) 3006
  73. SDouble Sqrt(const SDouble& x); // sqrt(x) 3007
  74. // cubic root
  75. SDouble Cbrt(const SDouble& x); // x^(1/3) 3010
  76. // reciporocal cubic root
  77. SDouble RecCbrt(const SDouble& x); // x^(-1/3) 3011
  78. // quartic root
  79. SDouble Qrrt(const SDouble& a); // x^(1/4) 3012
  80. // reciporocal quartic root
  81. SDouble RecQrrt(const SDouble& a); // x^(-1/4) 3013
  82. // exponential and hyperbolic functions
  83. SDouble Expower(long n); // n-th power of e 3302
  84. /**************************************************************************
  85. trigonometric functions calculated by DRAIDX
  86. **************************************************************************/
  87. SDouble Acos(const SDouble& x); // 3101
  88. SDouble Asin(const SDouble& x); // 3102
  89. inline SDouble Atan(const SDouble& x) { return AtanT(x); } // since version 2.21
  90. inline SDouble Cos(const SDouble& x) { return CosBS(x); } // since version 2.21
  91. inline SDouble Sin(const SDouble& x) { return SinBS(x); }
  92. inline SDouble Tan(const SDouble& x) { return TanBS(x); }
  93. // SDouble version of C's "double atan2(double y, double x)" arctan(y/x) 3110
  94. SDouble Atan2(const SDouble& y, const SDouble& x);
  95. /**************************************************************************
  96. It provides the exponential function exp(x).
  97. From the condition "exp(x) < DRADIX^DRADIX_EXP_MAX" the range of 'x' must be
  98. within
  99. |x| < DFIGURES*log(10)*DRADIX_EXP_MAX = 301759.17...
  100. to inline since version 2.21
  101. ***************************************************************************/
  102. inline SDouble Exp(const SDouble& x) {
  103. return ExpBS(x);
  104. }
  105. inline SDouble Cosh(const SDouble& x) {
  106. return CoshBS(x);
  107. }
  108. inline SDouble Sinh(const SDouble& x) {
  109. return SinhBS(x);
  110. }
  111. inline SDouble Tanh(const SDouble& x){
  112. return TanhBS(x);
  113. }
  114. // inverse hyperbolic functions
  115. SDouble Acosh(const SDouble& x); // 3309
  116. SDouble Asinh(const SDouble& x); // 3310
  117. SDouble Atanh(const SDouble& x); // 3311
  118. // logalithm functions
  119. inline SDouble Log(const SDouble& x) { return LogE(x); } // log(x)
  120. inline SDouble Log10(const SDouble& x){ return Log(x)*RecLog10(); } // log_10(x)
  121. /***********************************************************************
  122. Prototype declaration of constant functions has a formula
  123. SDouble func(const SDouble* userV = NULL,SDouble (*pfCalcFunc)() = fname);
  124. userV : a value given by user
  125. pfCalcFunc : default function, NULL is ok. A funcion which is faster than that
  126. of SN library can be used.
  127. See also "SetConstByFile()" in below.
  128. ************************************************************************/
  129. SDouble E(const SDouble* e = NULL,
  130. SDouble (*pfCalcFunc)(void) = SNE); // base of natural logarithm 3521
  131. SDouble Log10(const SDouble* ln10 = NULL,
  132. SDouble (*pfCalcFunc)(void) = SNLog10); // log(10.0) 3522
  133. SDouble Pi(const SDouble* pi = NULL,
  134. SDouble (*pfCalcFunc)(void) = SNPi); // pi 3523
  135. SDouble M2Pi(void); // 2/pi 3506
  136. SDouble RecE(void); // 1/e 3507
  137. SDouble MPi2(void); // pi/2 3511
  138. SDouble MPi4(void); // pi/4 3550
  139. SDouble RecLog10(void); // 1/log(10.0) 3551
  140. /////////////////// SComplex class /////////////////
  141. #ifndef SCOMPLEX_H
  142. #include "scomplex.h"
  143. #endif // SCOMPLEX_H
  144. inline SComplex Conj(const SComplex& z) { return z.Conj(); } // to inline since version 2.21
  145. inline SDouble Cabs(const SComplex& z) { return Sqrt(z.Norm()); } // |z| most simple
  146. inline SDouble Arg(const SComplex& z) { return Atan2(z.Imag(), z.Real()); } // argument in xy plane
  147. inline SDouble Real(const SComplex& z){ return z.Real();} // real part
  148. inline SDouble Imag(const SComplex& z){ return z.Imag();} // imaginary part
  149. inline SDouble Re(const SComplex& z) { return z.Real();} // Re(z) real part
  150. inline SDouble Im(const SComplex& z) { return z.Imag();} // Im(z) imaginary part
  151. inline SDouble Norm(const SComplex& z){ return z.Norm(); } // square of the magnitude
  152. /*********************************
  153. SComplex mathematical functions
  154. Acooding to GCC's "complex.h" the initial letter
  155. of funtion name is 'C';
  156. **********************************/
  157. SComplex Csqrt(const SComplex& z); // sqrt(z) 9101
  158. SComplex Cpolar(const SDouble& r, const SDouble& angle); // r*{cos(angle)+i*sin(angle)}; 9102
  159. SComplex Cpow(const SComplex& z, int n); // z^n 9103
  160. SComplex Cpow(const SComplex& x, const SDouble& y);// x^y 9104
  161. SComplex Cpow(const SDouble& x, const SComplex& y);// x^y 9105
  162. SComplex Cpow(const SComplex& x, const SComplex& y);// x^y 9106
  163. SComplex Cexp(const SComplex& z); // exp z 9107
  164. SComplex Clog(const SComplex& z); // log z 9108
  165. SComplex Ccosh(const SComplex& z); // cosh z 9109
  166. SComplex Csinh(const SComplex& z); // sinh z 9110
  167. SComplex Ctanh(const SComplex& z); // tanh z 9111
  168. SComplex Ccos(const SComplex& z); // cos z 9112
  169. SComplex Csin(const SComplex& z); // sin z 9113
  170. SComplex Ctan(const SComplex& z); // tan z 9114
  171. //The following are added since ver. 2.18.
  172. SComplex Cacos(const SComplex& z); // arccos z 9115
  173. SComplex Casin(const SComplex& z); // arcsin z 9116
  174. SComplex Catan(const SComplex& z); // arctan z 9117
  175. SComplex Cacosh(const SComplex& z); // arcosh z 9118
  176. SComplex Casinh(const SComplex& z); // arsinh z 9119
  177. SComplex Catanh(const SComplex& z); // artanh z 9120
  178. // See a sample program "smpcmpx3.cpp".
  179. // complex <double> complexC(const SComplex& z); //9151 ver. 2.17
  180. // conversion SComplex to complex <double>
  181. inline complex <double> complexC(const SComplex& z){
  182. return complex <double>(doubleD(z.Real()), doubleD(z.Imag()));
  183. }
  184. #endif // SN_MATH_H