/Python/pymath.c

http://unladen-swallow.googlecode.com/ · C · 247 lines · 150 code · 18 blank · 79 comment · 38 complexity · d92ee247d5047f14cd57bf2169354d21 MD5 · raw file

  1. #include "Python.h"
  2. #ifdef X87_DOUBLE_ROUNDING
  3. /* On x86 platforms using an x87 FPU, this function is called from the
  4. Py_FORCE_DOUBLE macro (defined in pymath.h) to force a floating-point
  5. number out of an 80-bit x87 FPU register and into a 64-bit memory location,
  6. thus rounding from extended precision to double precision. */
  7. double _Py_force_double(double x)
  8. {
  9. volatile double y;
  10. y = x;
  11. return y;
  12. }
  13. #endif
  14. #ifndef HAVE_HYPOT
  15. double hypot(double x, double y)
  16. {
  17. double yx;
  18. x = fabs(x);
  19. y = fabs(y);
  20. if (x < y) {
  21. double temp = x;
  22. x = y;
  23. y = temp;
  24. }
  25. if (x == 0.)
  26. return 0.;
  27. else {
  28. yx = y/x;
  29. return x*sqrt(1.+yx*yx);
  30. }
  31. }
  32. #endif /* HAVE_HYPOT */
  33. #ifndef HAVE_COPYSIGN
  34. double
  35. copysign(double x, double y)
  36. {
  37. /* use atan2 to distinguish -0. from 0. */
  38. if (y > 0. || (y == 0. && atan2(y, -1.) > 0.)) {
  39. return fabs(x);
  40. } else {
  41. return -fabs(x);
  42. }
  43. }
  44. #endif /* HAVE_COPYSIGN */
  45. #ifndef HAVE_LOG1P
  46. #include <float.h>
  47. double
  48. log1p(double x)
  49. {
  50. /* For x small, we use the following approach. Let y be the nearest
  51. float to 1+x, then
  52. 1+x = y * (1 - (y-1-x)/y)
  53. so log(1+x) = log(y) + log(1-(y-1-x)/y). Since (y-1-x)/y is tiny,
  54. the second term is well approximated by (y-1-x)/y. If abs(x) >=
  55. DBL_EPSILON/2 or the rounding-mode is some form of round-to-nearest
  56. then y-1-x will be exactly representable, and is computed exactly
  57. by (y-1)-x.
  58. If abs(x) < DBL_EPSILON/2 and the rounding mode is not known to be
  59. round-to-nearest then this method is slightly dangerous: 1+x could
  60. be rounded up to 1+DBL_EPSILON instead of down to 1, and in that
  61. case y-1-x will not be exactly representable any more and the
  62. result can be off by many ulps. But this is easily fixed: for a
  63. floating-point number |x| < DBL_EPSILON/2., the closest
  64. floating-point number to log(1+x) is exactly x.
  65. */
  66. double y;
  67. if (fabs(x) < DBL_EPSILON/2.) {
  68. return x;
  69. } else if (-0.5 <= x && x <= 1.) {
  70. /* WARNING: it's possible than an overeager compiler
  71. will incorrectly optimize the following two lines
  72. to the equivalent of "return log(1.+x)". If this
  73. happens, then results from log1p will be inaccurate
  74. for small x. */
  75. y = 1.+x;
  76. return log(y)-((y-1.)-x)/y;
  77. } else {
  78. /* NaNs and infinities should end up here */
  79. return log(1.+x);
  80. }
  81. }
  82. #endif /* HAVE_LOG1P */
  83. /*
  84. * ====================================================
  85. * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  86. *
  87. * Developed at SunPro, a Sun Microsystems, Inc. business.
  88. * Permission to use, copy, modify, and distribute this
  89. * software is freely granted, provided that this notice
  90. * is preserved.
  91. * ====================================================
  92. */
  93. static const double ln2 = 6.93147180559945286227E-01;
  94. static const double two_pow_m28 = 3.7252902984619141E-09; /* 2**-28 */
  95. static const double two_pow_p28 = 268435456.0; /* 2**28 */
  96. static const double zero = 0.0;
  97. /* asinh(x)
  98. * Method :
  99. * Based on
  100. * asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ]
  101. * we have
  102. * asinh(x) := x if 1+x*x=1,
  103. * := sign(x)*(log(x)+ln2)) for large |x|, else
  104. * := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else
  105. * := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2)))
  106. */
  107. #ifndef HAVE_ASINH
  108. double
  109. asinh(double x)
  110. {
  111. double w;
  112. double absx = fabs(x);
  113. if (Py_IS_NAN(x) || Py_IS_INFINITY(x)) {
  114. return x+x;
  115. }
  116. if (absx < two_pow_m28) { /* |x| < 2**-28 */
  117. return x; /* return x inexact except 0 */
  118. }
  119. if (absx > two_pow_p28) { /* |x| > 2**28 */
  120. w = log(absx)+ln2;
  121. }
  122. else if (absx > 2.0) { /* 2 < |x| < 2**28 */
  123. w = log(2.0*absx + 1.0 / (sqrt(x*x + 1.0) + absx));
  124. }
  125. else { /* 2**-28 <= |x| < 2= */
  126. double t = x*x;
  127. w = log1p(absx + t / (1.0 + sqrt(1.0 + t)));
  128. }
  129. return copysign(w, x);
  130. }
  131. #endif /* HAVE_ASINH */
  132. /* acosh(x)
  133. * Method :
  134. * Based on
  135. * acosh(x) = log [ x + sqrt(x*x-1) ]
  136. * we have
  137. * acosh(x) := log(x)+ln2, if x is large; else
  138. * acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else
  139. * acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1.
  140. *
  141. * Special cases:
  142. * acosh(x) is NaN with signal if x<1.
  143. * acosh(NaN) is NaN without signal.
  144. */
  145. #ifndef HAVE_ACOSH
  146. double
  147. acosh(double x)
  148. {
  149. if (Py_IS_NAN(x)) {
  150. return x+x;
  151. }
  152. if (x < 1.) { /* x < 1; return a signaling NaN */
  153. errno = EDOM;
  154. #ifdef Py_NAN
  155. return Py_NAN;
  156. #else
  157. return (x-x)/(x-x);
  158. #endif
  159. }
  160. else if (x >= two_pow_p28) { /* x > 2**28 */
  161. if (Py_IS_INFINITY(x)) {
  162. return x+x;
  163. } else {
  164. return log(x)+ln2; /* acosh(huge)=log(2x) */
  165. }
  166. }
  167. else if (x == 1.) {
  168. return 0.0; /* acosh(1) = 0 */
  169. }
  170. else if (x > 2.) { /* 2 < x < 2**28 */
  171. double t = x*x;
  172. return log(2.0*x - 1.0 / (x + sqrt(t - 1.0)));
  173. }
  174. else { /* 1 < x <= 2 */
  175. double t = x - 1.0;
  176. return log1p(t + sqrt(2.0*t + t*t));
  177. }
  178. }
  179. #endif /* HAVE_ACOSH */
  180. /* atanh(x)
  181. * Method :
  182. * 1.Reduced x to positive by atanh(-x) = -atanh(x)
  183. * 2.For x>=0.5
  184. * 1 2x x
  185. * atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
  186. * 2 1 - x 1 - x
  187. *
  188. * For x<0.5
  189. * atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
  190. *
  191. * Special cases:
  192. * atanh(x) is NaN if |x| >= 1 with signal;
  193. * atanh(NaN) is that NaN with no signal;
  194. *
  195. */
  196. #ifndef HAVE_ATANH
  197. double
  198. atanh(double x)
  199. {
  200. double absx;
  201. double t;
  202. if (Py_IS_NAN(x)) {
  203. return x+x;
  204. }
  205. absx = fabs(x);
  206. if (absx >= 1.) { /* |x| >= 1 */
  207. errno = EDOM;
  208. #ifdef Py_NAN
  209. return Py_NAN;
  210. #else
  211. return x/zero;
  212. #endif
  213. }
  214. if (absx < two_pow_m28) { /* |x| < 2**-28 */
  215. return x;
  216. }
  217. if (absx < 0.5) { /* |x| < 0.5 */
  218. t = absx+absx;
  219. t = 0.5 * log1p(t + t*absx / (1.0 - absx));
  220. }
  221. else { /* 0.5 <= |x| <= 1.0 */
  222. t = 0.5 * log1p((absx + absx) / (1.0 - absx));
  223. }
  224. return copysign(t, x);
  225. }
  226. #endif /* HAVE_ATANH */