Mark Dickinson added the comment:
Here is (quite a large) patch, cmath.patch, that fixes a variety of
problems in the current cmath module. A summary of the changes:
* sqrt, log, acos, acosh, asin, asinh, atan, atanh no longer produce
overflow errors for very large inputs
* exp, cos, sin, tan, cosh, sinh, tanh produce valid answers in some cases
where they incorrectly overflowed before.
* sqrt and log are more accurate for tiny numbers
* numeric problems in some functions (especially asinh and asin) should
have been fixed
* the branch cuts for asinh have been moved to the standard positions
* the direction of continuity for the branch cuts of tan, tanh have been
fixed
* on systems with full hardware and system support for signed zeros (most
modern systems), functions with a branch cut are continuous on both sides
of the branch cut. (As recommended by the C99 standard, amongst others.)
The patch includes documentation updates, but there are still no tests.
(My current tests make heavy use of the MPFR library, and assume IEEE-754
format doubles, so need to be seriously reworked.) The tests are on my to-
do list, but I'm unlikely to find time for them before the new year. In
the meantime, I'd very much appreciate any feedback on this patch, if
anyone has the time/energy/inclination to look at it. (Andreas: are you
in a position to give this a test-run?)
Added file: http://bugs.python.org/file8807/cmath.patch
__________________________________
Tracker <[EMAIL PROTECTED]>
<http://bugs.python.org/issue1381>
__________________________________
Index: configure
===================================================================
--- configure (revision 59184)
+++ configure (working copy)
@@ -1,5 +1,5 @@
#! /bin/sh
-# From configure.in Revision: 58653 .
+# From configure.in Revision: 58784 .
# Guess values for system-dependent variables and create Makefiles.
# Generated by GNU Autoconf 2.61 for python 2.6.
#
@@ -15226,11 +15226,14 @@
-for ac_func in alarm bind_textdomain_codeset chflags chown clock confstr \
- ctermid execv fork fpathconf ftime ftruncate \
+
+
+
+for ac_func in alarm asinh bind_textdomain_codeset chflags chown clock \
+ confstr copysign ctermid execv fork fpathconf ftime ftruncate \
gai_strerror getgroups getlogin getloadavg getpeername getpgid getpid \
getpriority getpwent getspnam getspent getsid getwd \
- kill killpg lchflags lchown lstat mkfifo mknod mktime \
+ kill killpg lchflags lchown log1p lstat mkfifo mknod mktime \
mremap nice pathconf pause plock poll pthread_init \
putenv readlink realpath \
select setegid seteuid setgid \
Index: configure.in
===================================================================
--- configure.in (revision 59184)
+++ configure.in (working copy)
@@ -2303,11 +2303,11 @@
AC_MSG_RESULT(MACHDEP_OBJS)
# checks for library functions
-AC_CHECK_FUNCS(alarm bind_textdomain_codeset chflags chown clock confstr \
- ctermid execv fork fpathconf ftime ftruncate \
+AC_CHECK_FUNCS(alarm asinh bind_textdomain_codeset chflags chown clock \
+ confstr copysign ctermid execv fork fpathconf ftime ftruncate \
gai_strerror getgroups getlogin getloadavg getpeername getpgid getpid \
getpriority getpwent getspnam getspent getsid getwd \
- kill killpg lchflags lchown lstat mkfifo mknod mktime \
+ kill killpg lchflags lchown log1p lstat mkfifo mknod mktime \
mremap nice pathconf pause plock poll pthread_init \
putenv readlink realpath \
select setegid seteuid setgid \
Index: Doc/library/cmath.rst
===================================================================
--- Doc/library/cmath.rst (revision 59184)
+++ Doc/library/cmath.rst (working copy)
@@ -14,6 +14,15 @@
floating-point number, respectively, and the function is then applied to the
result of the conversion.
+.. note::
+
+ On platforms with hardware and system-level support for signed
+ zeros, functions involving branch cuts are continuous on *both*
+ sides of the branch cut: the sign of the zero distinguishes one
+ side of the branch cut from the other. On platforms that do not
+ support signed zeros the continuity is as specified below.
+
+
The functions are:
@@ -37,32 +46,37 @@
.. function:: asinh(x)
- Return the hyperbolic arc sine of *x*. There are two branch cuts, extending
- left from ``±1j`` to ``±âj``, both continuous from above. These branch
cuts
- should be considered a bug to be corrected in a future release. The correct
- branch cuts should extend along the imaginary axis, one from ``1j`` up to
- ``âj`` and continuous from the right, and one from ``-1j`` down to
``-âj``
- and continuous from the left.
+ Return the hyperbolic arc sine of *x*. There are two branch cuts:
+ One extends from ``1j`` along the imaginary axis to ``âj``,
+ continuous from the right. The other extends from ``-1j`` along
+ the imaginary axis to ``-âj``, continuous from the left.
+ .. versionchanged:: 2.6
+ branch cuts moved to match those recommended by the C99 standard
+
.. function:: atan(x)
Return the arc tangent of *x*. There are two branch cuts: One extends from
- ``1j`` along the imaginary axis to ``âj``, continuous from the left. The
+ ``1j`` along the imaginary axis to ``âj``, continuous from the right. The
other extends from ``-1j`` along the imaginary axis to ``-âj``, continuous
- from the left. (This should probably be changed so the upper cut becomes
- continuous from the other side.)
+ from the left.
+ .. versionchanged:: 2.6
+ direction of continuity of upper cut reversed
+
.. function:: atanh(x)
Return the hyperbolic arc tangent of *x*. There are two branch cuts: One
- extends from ``1`` along the real axis to ``â``, continuous from above.
The
+ extends from ``1`` along the real axis to ``â``, continuous from below.
The
other extends from ``-1`` along the real axis to ``-â``, continuous from
- above. (This should probably be changed so the right cut becomes continuous
- from the other side.)
+ above.
+ .. versionchanged:: 2.6
+ direction of continuity of right cut reversed
+
.. function:: cos(x)
Return the cosine of *x*.
Index: Modules/cmathmodule.c
===================================================================
--- Modules/cmathmodule.c (revision 59184)
+++ Modules/cmathmodule.c (working copy)
@@ -5,29 +5,172 @@
#include "Python.h"
#ifndef M_PI
-#define M_PI (3.141592653589793239)
+#define M_PI (3.141592653589793238)
#endif
-/* First, the C functions that do the real work */
+#ifndef M_PI_2
+#define M_PI_2 (1.570796326794896619)
+#endif
-/* constants */
-static Py_complex c_one = {1., 0.};
-static Py_complex c_half = {0.5, 0.};
-static Py_complex c_i = {0., 1.};
-static Py_complex c_halfi = {0., 0.5};
+#ifndef M_E
+#define M_E (2.718281828459045235)
+#endif
+#ifndef M_LN2
+#define M_LN2 (0.6931471805599453094) /* natural log of 2 */
+#endif
+
+#ifndef M_LN4
+#define M_LN4 (1.386294361119890619) /* natural log of 4 */
+#endif
+
+#ifndef M_LN10
+#define M_LN10 (2.302585092994045684) /* natural log of 10 */
+#endif
+
+
+/* Machine-dependent constants used as cutoff points for some of the
+ algorithms below. The values chosen below should work for DEC VAX, IEEE
+ and IBM S/390 double-precision formats.
+
+ LARGE_DOUBLE is a large positive number used to avoid spurious overflow in
+ the sqrt, log, inverse trig and inverse hyperbolic trig functions. It
should be large
+ enough that 1/sqrt(LARGE_DOUBLE) is significantly smaller than the machine
+ epsilon (around 1e-16 on IEEE machines), and small enough that
+ 4*LARGE_DOUBLE does not overflow. SQRT_LARGE_DOUBLE is (approximately) the
+ square root of LARGE_DOUBLE; R_LARGE_DOUBLE is (approximately) the
+ reciprocal of LARGE_DOUBLE.
+
+ EXP_CUTOFF is used in the evaluation of exp, cos, cosh, sin, sinh, tan,
+ tanh to avoid unecessary overflow. It should be small enough that
+ exp(EXP_CUTOFF) is representable.
+
+ LOG1P_TAYLOR_CUTOFF is used to decide when to use a Taylor-series expansion
+ for the log1p function. It should be chosen so that LOG1P_TAYLOR_CUTOFF is
+ larger than the machine epsilon, but its square is smaller than the machine
+ epsilon.
+ */
+
+#define LARGE_DOUBLE (1e300)
+#define R_LARGE_DOUBLE (1./LARGE_DOUBLE)
+#define SQRT_LARGE_DOUBLE (1e150)
+#define R_SQRT_LARGE_DOUBLE (1./SQRT_LARGE_DOUBLE)
+#define EXP_CUTOFF (700.)
+#define LOG1P_TAYLOR_CUTOFF (1e-10)
+
+
/* forward declarations */
-static Py_complex c_log(Py_complex);
-static Py_complex c_prodi(Py_complex);
+static Py_complex c_asinh(Py_complex);
+static Py_complex c_atanh(Py_complex);
+static Py_complex c_cosh(Py_complex);
+static Py_complex c_sinh(Py_complex);
static Py_complex c_sqrt(Py_complex);
+static Py_complex c_tanh(Py_complex);
static PyObject * math_error(void);
+/* the copysign, log1p and asinh functions are all part of the C99 standard,
+ but not part of the C89 standard. Here are substitutes for these three
+ functions, to be used if the math library doesn't already have them. */
+#ifndef HAVE_COPYSIGN
+static double
+copysign(double x, double y)
+{
+ /* use atan2 to distinguish -0. from 0. */
+ if (y > 0. || (y == 0. && atan2(y, -1.) > 0.)) {
+ return fabs(x);
+ } else {
+ return -fabs(x);
+ }
+}
+#endif
+
+#ifndef HAVE_LOG1P
+static double
+log1p(double x)
+{
+ double y;
+ if (fabs(x) < LOG1P_TAYLOR_CUTOFF) {
+ /* for tiny arguments, use two terms of the Taylor series. */
+ return x*(1.-x/2.);
+ } else if (-0.5 <= x && x <= 1.) {
+ /* WARNING: it's possible than an errant compiler will
+ incorrectly optimize the following two lines to the
+ equivalent of "return log(1.+x)". If this happens, then
+ results from log1p will be inaccurate for small x. */
+ y = 1.+x;
+ return log(y)-((y-1.)-x)/y;
+ } else {
+ /* NaNs and infinities should end up here */
+ return log(1.+x);
+ }
+}
+#endif
+
+#ifndef HAVE_ASINH
+static double
+asinh(double x)
+{
+ /* asinh is an odd function, so we can use
+
+ asinh(x) = copysign(asinh(abs(x)), x)
+
+ to reduce to the case when x is nonnegative. Here, the usual
+ formula is:
+
+ asinh(x) = log(x + hypot(1, x))
+
+ For x small, the argument of the log is close to 1, resulting in
+ a loss of precision; in this case we use the formula
+
+ asinh(x) = log1p(x + x*x/(hypot(1, x)+1))
+
+ For x very large there's a danger of overflow in computing x +
hypot(1, x).
+ In this case we use the approximation
+
+ asinh(x) ~ log(2*x) = log(x) + log(2).
+
+ */
+
+ double ax;
+ ax = fabs(x);
+ if (ax > LARGE_DOUBLE) {
+ return copysign(log(ax) + M_LN2, x);
+ } else if (ax < 1.) {
+ return copysign(log1p(ax*(1.+ax/(hypot(1., ax)+1.))), x);
+ } else {
+ return copysign(log(ax + hypot(1., ax)), x);
+ }
+}
+#endif
+
+/* First, the C functions that do the real work */
+
static Py_complex
-c_acos(Py_complex x)
+c_acos(Py_complex z)
{
- return c_neg(c_prodi(c_log(c_sum(x,c_prod(c_i,
- c_sqrt(c_diff(c_one,c_prod(x,x))))))));
+ Py_complex s1, s2, r;
+ if (fabs(z.real) > LARGE_DOUBLE || fabs(z.imag) > LARGE_DOUBLE) {
+ /* avoid unnecessary overflow for large arguments */
+ r.real = atan2(fabs(z.imag), z.real);
+ /* split into cases to make sure that the branch cut has the
+ correct continuity on systems with unsigned zeros */
+ if (z.real < 0.) {
+ r.imag = -copysign(log(hypot(z.real/2., z.imag/2.)) +
M_LN4, z.imag);
+ } else {
+ r.imag = copysign(log(hypot(z.real/2., z.imag/2.)) +
M_LN4, -z.imag);
+ }
+ } else {
+ s1.real = 1.-z.real;
+ s1.imag = -z.imag;
+ s1 = c_sqrt(s1);
+ s2.real = 1.+z.real;
+ s2.imag = z.imag;
+ s2 = c_sqrt(s2);
+ r.real = 2.*atan2(s1.real, s2.real);
+ r.imag = asinh(s2.real*s1.imag - s2.imag*s1.real);
+ }
+ return r;
}
PyDoc_STRVAR(c_acos_doc,
@@ -37,13 +180,25 @@
static Py_complex
-c_acosh(Py_complex x)
+c_acosh(Py_complex z)
{
- Py_complex z;
- z = c_sqrt(c_half);
- z = c_log(c_prod(z, c_sum(c_sqrt(c_sum(x,c_one)),
- c_sqrt(c_diff(x,c_one)))));
- return c_sum(z, z);
+ Py_complex s1, s2, r;
+
+ if (fabs(z.real) > LARGE_DOUBLE || fabs(z.imag) > LARGE_DOUBLE) {
+ /* avoid unnecessary overflow for large arguments */
+ r.real = log(hypot(z.real/2., z.imag/2.)) + M_LN4;
+ r.imag = atan2(z.imag, z.real);
+ } else {
+ s1.real = z.real - 1.;
+ s1.imag = z.imag;
+ s1 = c_sqrt(s1);
+ s2.real = z.real + 1.;
+ s2.imag = z.imag;
+ s2 = c_sqrt(s2);
+ r.real = asinh(s1.real*s2.real + s1.imag*s2.imag);
+ r.imag = 2.*atan2(s1.imag, s2.real);
+ }
+ return r;
}
PyDoc_STRVAR(c_acosh_doc,
@@ -53,14 +208,16 @@
static Py_complex
-c_asin(Py_complex x)
+c_asin(Py_complex z)
{
- /* -i * log[(sqrt(1-x**2) + i*x] */
- const Py_complex squared = c_prod(x, x);
- const Py_complex sqrt_1_minus_x_sq = c_sqrt(c_diff(c_one, squared));
- return c_neg(c_prodi(c_log(
- c_sum(sqrt_1_minus_x_sq, c_prodi(x))
- ) ) );
+ /* asin(z) = -i asinh(iz) */
+ Py_complex s, r;
+ s.real = -z.imag;
+ s.imag = z.real;
+ s = c_asinh(s);
+ r.real = s.imag;
+ r.imag = -s.real;
+ return r;
}
PyDoc_STRVAR(c_asin_doc,
@@ -70,13 +227,28 @@
static Py_complex
-c_asinh(Py_complex x)
+c_asinh(Py_complex z)
{
- Py_complex z;
- z = c_sqrt(c_half);
- z = c_log(c_prod(z, c_sum(c_sqrt(c_sum(x, c_i)),
- c_sqrt(c_diff(x, c_i)))));
- return c_sum(z, z);
+ Py_complex s1, s2, r;
+
+ if (fabs(z.real) > LARGE_DOUBLE || fabs(z.imag) > LARGE_DOUBLE) {
+ if (z.imag >= 0.) {
+ r.real = copysign(log(hypot(z.real/2., z.imag/2.)) +
M_LN4, z.real);
+ } else {
+ r.real = -copysign(log(hypot(z.real/2., z.imag/2.)) +
M_LN4, -z.real);
+ }
+ r.imag = atan2(z.imag, fabs(z.real));
+ } else {
+ s1.real = 1.+z.imag;
+ s1.imag = -z.real;
+ s1 = c_sqrt(s1);
+ s2.real = 1.-z.imag;
+ s2.imag = z.real;
+ s2 = c_sqrt(s2);
+ r.real = asinh(s1.real*s2.imag-s2.real*s1.imag);
+ r.imag = atan2(z.imag, s1.real*s2.real-s1.imag*s2.imag);
+ }
+ return r;
}
PyDoc_STRVAR(c_asinh_doc,
@@ -86,9 +258,16 @@
static Py_complex
-c_atan(Py_complex x)
+c_atan(Py_complex z)
{
- return c_prod(c_halfi,c_log(c_quot(c_sum(c_i,x),c_diff(c_i,x))));
+ /* atan(z) = -i atanh(iz) */
+ Py_complex s, r;
+ s.real = -z.imag;
+ s.imag = z.real;
+ s = c_atanh(s);
+ r.real = s.imag;
+ r.imag = -s.real;
+ return r;
}
PyDoc_STRVAR(c_atan_doc,
@@ -98,9 +277,43 @@
static Py_complex
-c_atanh(Py_complex x)
+c_atanh(Py_complex z)
{
- return c_prod(c_half,c_log(c_quot(c_sum(c_one,x),c_diff(c_one,x))));
+ Py_complex r;
+ double ay, ratio;
+
+ /* use oddness to reduce to case where real part is >= 0. */
+ if (z.real < 0.) {
+ return c_neg(c_atanh(c_neg(z)));
+ }
+
+ ay = fabs(z.imag);
+ if (z.real > SQRT_LARGE_DOUBLE || ay > SQRT_LARGE_DOUBLE) {
+ if (ay < z.real) {
+ ratio = ay/z.real;
+ r.real = 1./(1.+ratio*ratio)/z.real;
+ } else {
+ ratio = z.real/ay;
+ r.real = ratio/(1.+ratio*ratio)/ay;
+ }
+ /* the two negations in the next line cancel each other out
+ except when working with unsigned zeros: they're there to
+ ensure that the branch cut has the correct continuity on
+ systems that don't support signed zeros */
+ r.imag = -copysign(M_PI_2, -z.imag);
+ } else if (z.real == 1. && ay < R_SQRT_LARGE_DOUBLE) {
+ /* C99 standard says: atanh(1+/-0.) should be inf +/- 0i */
+ r.real = -log(sqrt(ay)/sqrt(hypot(ay, 2.)));
+ if (ay == 0.) {
+ r.imag = z.imag;
+ } else {
+ r.imag = copysign(atan2(2., -ay)/2, z.imag);
+ }
+ } else {
+ r.real = log1p(4.*z.real/((1-z.real)*(1-z.real) + ay*ay))/4.;
+ r.imag = -atan2(-2.*z.imag, (1-z.real)*(1+z.real) - ay*ay)/2.;
+ }
+ return r;
}
PyDoc_STRVAR(c_atanh_doc,
@@ -110,11 +323,13 @@
static Py_complex
-c_cos(Py_complex x)
+c_cos(Py_complex z)
{
+ /* cos(z) = cosh(iz) */
Py_complex r;
- r.real = cos(x.real)*cosh(x.imag);
- r.imag = -sin(x.real)*sinh(x.imag);
+ r.real = -z.imag;
+ r.imag = z.real;
+ r = c_cosh(r);
return r;
}
@@ -125,11 +340,20 @@
static Py_complex
-c_cosh(Py_complex x)
+c_cosh(Py_complex z)
{
Py_complex r;
- r.real = cos(x.imag)*cosh(x.real);
- r.imag = sin(x.imag)*sinh(x.real);
+ double x_minus_one;
+
+ if (fabs(z.real) > EXP_CUTOFF) {
+ /* avoid corner cases involving unnecessary overflow */
+ x_minus_one = z.real - copysign(1., z.real);
+ r.real = cos(z.imag) * cosh(x_minus_one) * M_E;
+ r.imag = sin(z.imag) * sinh(x_minus_one) * M_E;
+ } else {
+ r.real = cos(z.imag) * cosh(z.real);
+ r.imag = sin(z.imag) * sinh(z.real);
+ }
return r;
}
@@ -140,12 +364,20 @@
static Py_complex
-c_exp(Py_complex x)
+c_exp(Py_complex z)
{
Py_complex r;
- double l = exp(x.real);
- r.real = l*cos(x.imag);
- r.imag = l*sin(x.imag);
+ double l;
+
+ if (z.real > EXP_CUTOFF) {
+ l = exp(z.real-1.);
+ r.real = l*cos(z.imag)*M_E;
+ r.imag = l*sin(z.imag)*M_E;
+ } else {
+ l = exp(z.real);
+ r.real = l*cos(z.imag);
+ r.imag = l*sin(z.imag);
+ }
return r;
}
@@ -156,23 +388,75 @@
static Py_complex
-c_log(Py_complex x)
+c_log(Py_complex z)
{
+ /*
+ The usual formula for the real part is log(hypot(z.real, z.imag)).
+ There are four situations where this formula is potentially
+ problematic:
+
+ (1) the absolute value of z is subnormal. Then hypot is subnormal,
+ so has fewer than the usual number of bits of accuracy, hence may
+ have large relative error. This then gives a large relative error
+ in the log. (But note that this is not a problem if either z.real
+ or z.imag is zero, in which case hypot(z.real, z.imag) is
+ represented exactly.)
+
+ (2) The absolute value of z overflows (e.g. when both z.real and
+ z.imag are within a factor of 1/sqrt(2) of the largest
+ representable double).
+
+ (3) the absolute value of z is close to 1. In this case it's
+ unreasonable to expect good accuracy: a change of 1ulp in the real
+ or imaginary part of z can result in a change of billions of ulps
+ in the correctly rounded answer. But there's one special case
+ where good accuracy *is* achievable with little effort, namely when
+ fabs(z.real) == 1 or fabs(z.imag) == 1.
+
+ (4) z = 0. The simplest thing to do here is to call the
+ floating-point log with an argument of 0, and let its behaviour
+ (returning -infinity, signaling a floating-point exception, setting
+ errno, or whatever) determine that of c_log. So the usual formula
+ is fine here.
+
+ All of situations (1)-(3) above can be resolved by using the
+ formula
+
+ log(hypot(x, y)) = log(y) + log1p((x/y)*(x/y))/2
+
+ for 0 < x <= y.
+
+ */
+
Py_complex r;
- double l = hypot(x.real,x.imag);
- r.imag = atan2(x.imag, x.real);
- r.real = log(l);
+ double ax, ay, am, ratio;
+
+ ax = fabs(z.real);
+ ay = fabs(z.imag);
+ if ((0. < ax && ax < R_LARGE_DOUBLE && 0. < ay && ay < R_LARGE_DOUBLE)
||
+ ax > LARGE_DOUBLE ||
+ ay > LARGE_DOUBLE ||
+ ax == 1. ||
+ ay == 1.) {
+ am = (ax > ay ? ax : ay); /* am = MAX(ax, ay) */
+ ratio = (ax > ay ? ay : ax)/am; /* ratio = MIN(ax, ay)/MAX(ax,
ay) */
+ r.real = log(am) + log1p(ratio*ratio)/2.;
+ } else {
+ r.real = log(hypot(ax, ay));
+ }
+ r.imag = atan2(z.imag, z.real);
return r;
}
static Py_complex
-c_log10(Py_complex x)
+c_log10(Py_complex z)
{
Py_complex r;
- double l = hypot(x.real,x.imag);
- r.imag = atan2(x.imag, x.real)/log(10.);
- r.real = log10(l);
+
+ r = c_log(z);
+ r.real = r.real / M_LN10;
+ r.imag = r.imag / M_LN10;
return r;
}
@@ -182,26 +466,19 @@
"Return the base-10 logarithm of x.");
-/* internal function not available from Python */
static Py_complex
-c_prodi(Py_complex x)
+c_sin(Py_complex z)
{
- Py_complex r;
- r.real = -x.imag;
- r.imag = x.real;
+ /* sin(z) = -i sin(iz) */
+ Py_complex s, r;
+ s.real = -z.imag;
+ s.imag = z.real;
+ s = c_sinh(s);
+ r.real = s.imag;
+ r.imag = -s.real;
return r;
}
-
-static Py_complex
-c_sin(Py_complex x)
-{
- Py_complex r;
- r.real = sin(x.real) * cosh(x.imag);
- r.imag = cos(x.real) * sinh(x.imag);
- return r;
-}
-
PyDoc_STRVAR(c_sin_doc,
"sin(x)\n"
"\n"
@@ -209,12 +486,21 @@
static Py_complex
-c_sinh(Py_complex x)
+c_sinh(Py_complex z)
{
Py_complex r;
- r.real = cos(x.imag) * sinh(x.real);
- r.imag = sin(x.imag) * cosh(x.real);
+ double x_minus_one;
+
+ if (fabs(z.real) > EXP_CUTOFF) {
+ x_minus_one = z.real - copysign(1., z.real);
+ r.real = cos(z.imag) * sinh(x_minus_one) * M_E;
+ r.imag = sin(z.imag) * cosh(x_minus_one) * M_E;
+ } else {
+ r.real = cos(z.imag) * sinh(z.real);
+ r.imag = sin(z.imag) * cosh(z.real);
+ }
return r;
+
}
PyDoc_STRVAR(c_sinh_doc,
@@ -224,28 +510,66 @@
static Py_complex
-c_sqrt(Py_complex x)
+c_sqrt(Py_complex z)
{
+ /*
+ Method: use symmetries to reduce to the case when x = z.real and y
+ = z.imag are nonnegative. Then the real part of the result is
+ given by
+
+ s = sqrt((x + hypot(x, y))/2)
+
+ and the imaginary part is
+
+ d = (y/2)/s
+
+ If either x or y is very large then there's a risk of overflow in
+ computation of the expression x + hypot(x, y). We can avoid this
+ by rewriting the formula for s as:
+
+ s = 2*sqrt(x/8 + hypot(x/8, y/8))
+
+ This costs us two extra multiplications, but avoids the overhead of
+ checking for x and y large.
+
+ If both x and y are subnormal then hypot(x, y) may also be
+ subnormal, so will lack full precision. We solve this by rescaling
+ x and y by a sufficiently large power of 2 to ensure that x and y
+ are normal.
+ */
+
+
Py_complex r;
double s,d;
- if (x.real == 0. && x.imag == 0.)
- r = x;
- else {
- s = sqrt(0.5*(fabs(x.real) + hypot(x.real,x.imag)));
- d = 0.5*x.imag/s;
- if (x.real > 0.) {
- r.real = s;
- r.imag = d;
- }
- else if (x.imag >= 0.) {
- r.real = d;
- r.imag = s;
- }
- else {
- r.real = -d;
- r.imag = -s;
- }
+ double ax, ay;
+
+ if (z.real == 0. && z.imag == 0.) {
+ r.real = 0.;
+ r.imag = z.imag;
+ return r;
}
+
+ ax = fabs(z.real);
+ ay = fabs(z.imag);
+
+ if (ax < R_LARGE_DOUBLE && ay < R_LARGE_DOUBLE) {
+ /* here we catch cases where hypot(ax, ay) is subnormal */
+ ax = ldexp(ax, 65);
+ s = ldexp(sqrt(ax + hypot(ax, ldexp(ay, 65))), -33);
+ d = ay/(2.*s);
+ } else {
+ ax /= 8.;
+ s = 2.*sqrt(ax + hypot(ax, ay/8.));
+ d = ay/(2.*s);
+ }
+
+ if (z.real >= 0.) {
+ r.real = s;
+ r.imag = copysign(d, z.imag);
+ } else {
+ r.real = d;
+ r.imag = copysign(s, z.imag);
+ }
return r;
}
@@ -256,23 +580,15 @@
static Py_complex
-c_tan(Py_complex x)
+c_tan(Py_complex z)
{
- Py_complex r;
- double sr,cr,shi,chi;
- double rs,is,rc,ic;
- double d;
- sr = sin(x.real);
- cr = cos(x.real);
- shi = sinh(x.imag);
- chi = cosh(x.imag);
- rs = sr * chi;
- is = cr * shi;
- rc = cr * chi;
- ic = -sr * shi;
- d = rc*rc + ic * ic;
- r.real = (rs*rc + is*ic) / d;
- r.imag = (is*rc - rs*ic) / d;
+ /* tan(z) = -i tanh(iz) */
+ Py_complex s, r;
+ s.real = -z.imag;
+ s.imag = z.real;
+ s = c_tanh(s);
+ r.real = s.imag;
+ r.imag = -s.real;
return r;
}
@@ -283,23 +599,35 @@
static Py_complex
-c_tanh(Py_complex x)
+c_tanh(Py_complex z)
{
+ /* Formula:
+
+ tanh(x+iy) = (tanh(x)(1+tan(y)^2) + i tan(y)(1-tanh(x))^2) /
+ (1+tan(y)^2 tanh(x)^2)
+
+ To avoid excessive roundoff error, 1-tanh(x)^2 is better computed
+ as 1/cosh(x)^2. When abs(x) is large, we approximate 1-tanh(x)^2
+ by 4 exp(-2*x) instead, to avoid possible overflow in the
+ computation of cosh(x).
+
+ */
+
Py_complex r;
- double si,ci,shr,chr;
- double rs,is,rc,ic;
- double d;
- si = sin(x.imag);
- ci = cos(x.imag);
- shr = sinh(x.real);
- chr = cosh(x.real);
- rs = ci * shr;
- is = si * chr;
- rc = ci * chr;
- ic = si * shr;
- d = rc*rc + ic*ic;
- r.real = (rs*rc + is*ic) / d;
- r.imag = (is*rc - rs*ic) / d;
+ double tx, ty, cx, txty, denom;
+
+ if (fabs(z.real) > EXP_CUTOFF) {
+ r.real = copysign(1., z.real);
+ r.imag = 2.*sin(2.*z.imag)*exp(-2.*fabs(z.real));
+ } else {
+ tx = tanh(z.real);
+ ty = tan(z.imag);
+ cx = 1./cosh(z.real);
+ txty = tx*ty;
+ denom = 1. + txty*txty;
+ r.real = tx*(1.+ty*ty)/denom;
+ r.imag = ((ty/denom)*cx)*cx;
+ }
return r;
}
@@ -323,9 +651,9 @@
if (PyTuple_GET_SIZE(args) == 2)
x = c_quot(x, c_log(y));
PyFPE_END_PROTECT(x)
+ Py_ADJUST_ERANGE2(x.real, x.imag);
if (errno != 0)
return math_error();
- Py_ADJUST_ERANGE2(x.real, x.imag);
return PyComplex_FromCComplex(x);
}
Index: pyconfig.h.in
===================================================================
--- pyconfig.h.in (revision 59184)
+++ pyconfig.h.in (working copy)
@@ -37,6 +37,9 @@
/* Define this if your time.h defines altzone. */
#undef HAVE_ALTZONE
+/* Define to 1 if you have the `asinh' function. */
+#undef HAVE_ASINH
+
/* Define to 1 if you have the <asm/types.h> header file. */
#undef HAVE_ASM_TYPES_H
@@ -85,6 +88,9 @@
/* Define to 1 if you have the <conio.h> header file. */
#undef HAVE_CONIO_H
+/* Define to 1 if you have the `copysign' function. */
+#undef HAVE_COPYSIGN
+
/* Define to 1 if you have the `ctermid' function. */
#undef HAVE_CTERMID
@@ -333,6 +339,9 @@
/* Define to 1 if you have the <linux/netlink.h> header file. */
#undef HAVE_LINUX_NETLINK_H
+/* Define to 1 if you have the `log1p' function. */
+#undef HAVE_LOG1P
+
/* Define this if you have the type long long. */
#undef HAVE_LONG_LONG
_______________________________________________
Python-bugs-list mailing list
Unsubscribe:
http://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com
