diff options
Diffstat (limited to 'Modules')
| -rw-r--r-- | Modules/cmathmodule.c | 1019 | ||||
| -rw-r--r-- | Modules/mathmodule.c | 449 |
2 files changed, 1235 insertions, 233 deletions
diff --git a/Modules/cmathmodule.c b/Modules/cmathmodule.c index ec48ce8..8e3c31e 100644 --- a/Modules/cmathmodule.c +++ b/Modules/cmathmodule.c @@ -3,31 +3,172 @@ /* much code borrowed from mathmodule.c */ #include "Python.h" +/* we need DBL_MAX, DBL_MIN, DBL_EPSILON, DBL_MANT_DIG and FLT_RADIX from + float.h. We assume that FLT_RADIX is either 2 or 16. */ +#include <float.h> -#ifndef M_PI -#define M_PI (3.141592653589793239) +#if (FLT_RADIX != 2 && FLT_RADIX != 16) +#error "Modules/cmathmodule.c expects FLT_RADIX to be 2 or 16" #endif -/* First, the C functions that do the real work */ +#ifndef M_LN2 +#define M_LN2 (0.6931471805599453094) /* natural log of 2 */ +#endif + +#ifndef M_LN10 +#define M_LN10 (2.302585092994045684) /* natural log of 10 */ +#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}; +/* + CM_LARGE_DOUBLE is used to avoid spurious overflow in the sqrt, log, + inverse trig and inverse hyperbolic trig functions. Its log is used in the + evaluation of exp, cos, cosh, sin, sinh, tan, and tanh to avoid unecessary + overflow. + */ + +#define CM_LARGE_DOUBLE (DBL_MAX/4.) +#define CM_SQRT_LARGE_DOUBLE (sqrt(CM_LARGE_DOUBLE)) +#define CM_LOG_LARGE_DOUBLE (log(CM_LARGE_DOUBLE)) +#define CM_SQRT_DBL_MIN (sqrt(DBL_MIN)) + +/* + CM_SCALE_UP is an odd integer chosen such that multiplication by + 2**CM_SCALE_UP is sufficient to turn a subnormal into a normal. + CM_SCALE_DOWN is (-(CM_SCALE_UP+1)/2). These scalings are used to compute + square roots accurately when the real and imaginary parts of the argument + are subnormal. +*/ + +#if FLT_RADIX==2 +#define CM_SCALE_UP (2*(DBL_MANT_DIG/2) + 1) +#elif FLT_RADIX==16 +#define CM_SCALE_UP (4*DBL_MANT_DIG+1) +#endif +#define CM_SCALE_DOWN (-(CM_SCALE_UP+1)/2) /* 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); +/* Code to deal with special values (infinities, NaNs, etc.). */ + +/* special_type takes a double and returns an integer code indicating + the type of the double as follows: +*/ + +enum special_types { + ST_NINF, /* 0, negative infinity */ + ST_NEG, /* 1, negative finite number (nonzero) */ + ST_NZERO, /* 2, -0. */ + ST_PZERO, /* 3, +0. */ + ST_POS, /* 4, positive finite number (nonzero) */ + ST_PINF, /* 5, positive infinity */ + ST_NAN, /* 6, Not a Number */ +}; + +static enum special_types +special_type(double d) +{ + if (Py_IS_FINITE(d)) { + if (d != 0) { + if (copysign(1., d) == 1.) + return ST_POS; + else + return ST_NEG; + } + else { + if (copysign(1., d) == 1.) + return ST_PZERO; + else + return ST_NZERO; + } + } + if (Py_IS_NAN(d)) + return ST_NAN; + if (copysign(1., d) == 1.) + return ST_PINF; + else + return ST_NINF; +} + +#define SPECIAL_VALUE(z, table) \ + if (!Py_IS_FINITE((z).real) || !Py_IS_FINITE((z).imag)) { \ + errno = 0; \ + return table[special_type((z).real)] \ + [special_type((z).imag)]; \ + } + +#define P Py_MATH_PI +#define P14 0.25*Py_MATH_PI +#define P12 0.5*Py_MATH_PI +#define P34 0.75*Py_MATH_PI +#ifdef MS_WINDOWS +/* On Windows HUGE_VAL is an extern variable and not a constant. Since the + special value arrays need a constant we have to roll our own infinity + and nan. */ +# define INF (DBL_MAX*DBL_MAX) +# define N (INF*0.) +#else +# define INF Py_HUGE_VAL +# define N Py_NAN +#endif /* MS_WINDOWS */ +#define U -9.5426319407711027e33 /* unlikely value, used as placeholder */ + +/* First, the C functions that do the real work. Each of the c_* + functions computes and returns the C99 Annex G recommended result + and also sets errno as follows: errno = 0 if no floating-point + exception is associated with the result; errno = EDOM if C99 Annex + G recommends raising divide-by-zero or invalid for this result; and + errno = ERANGE where the overflow floating-point signal should be + raised. +*/ + +static Py_complex acos_special_values[7][7] = { + {{P34,INF},{P,INF}, {P,INF}, {P,-INF}, {P,-INF}, {P34,-INF},{N,INF}}, + {{P12,INF},{U,U}, {U,U}, {U,U}, {U,U}, {P12,-INF},{N,N}}, + {{P12,INF},{U,U}, {P12,0.},{P12,-0.},{U,U}, {P12,-INF},{P12,N}}, + {{P12,INF},{U,U}, {P12,0.},{P12,-0.},{U,U}, {P12,-INF},{P12,N}}, + {{P12,INF},{U,U}, {U,U}, {U,U}, {U,U}, {P12,-INF},{N,N}}, + {{P14,INF},{0.,INF},{0.,INF},{0.,-INF},{0.,-INF},{P14,-INF},{N,INF}}, + {{N,INF}, {N,N}, {N,N}, {N,N}, {N,N}, {N,-INF}, {N,N}} +}; 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; + + SPECIAL_VALUE(z, acos_special_values); + + if (fabs(z.real) > CM_LARGE_DOUBLE || fabs(z.imag) > CM_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_LN2*2., z.imag); + } else { + r.imag = copysign(log(hypot(z.real/2., z.imag/2.)) + + M_LN2*2., -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); + } + errno = 0; + return r; } PyDoc_STRVAR(c_acos_doc, @@ -36,14 +177,39 @@ PyDoc_STRVAR(c_acos_doc, "Return the arc cosine of x."); +static Py_complex acosh_special_values[7][7] = { + {{INF,-P34},{INF,-P}, {INF,-P}, {INF,P}, {INF,P}, {INF,P34},{INF,N}}, + {{INF,-P12},{U,U}, {U,U}, {U,U}, {U,U}, {INF,P12},{N,N}}, + {{INF,-P12},{U,U}, {0.,-P12},{0.,P12},{U,U}, {INF,P12},{N,N}}, + {{INF,-P12},{U,U}, {0.,-P12},{0.,P12},{U,U}, {INF,P12},{N,N}}, + {{INF,-P12},{U,U}, {U,U}, {U,U}, {U,U}, {INF,P12},{N,N}}, + {{INF,-P14},{INF,-0.},{INF,-0.},{INF,0.},{INF,0.},{INF,P14},{INF,N}}, + {{INF,N}, {N,N}, {N,N}, {N,N}, {N,N}, {INF,N}, {N,N}} +}; + 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; + + SPECIAL_VALUE(z, acosh_special_values); + + if (fabs(z.real) > CM_LARGE_DOUBLE || fabs(z.imag) > CM_LARGE_DOUBLE) { + /* avoid unnecessary overflow for large arguments */ + r.real = log(hypot(z.real/2., z.imag/2.)) + M_LN2*2.; + 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); + } + errno = 0; + return r; } PyDoc_STRVAR(c_acosh_doc, @@ -53,14 +219,16 @@ PyDoc_STRVAR(c_acosh_doc, 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, @@ -69,14 +237,44 @@ PyDoc_STRVAR(c_asin_doc, "Return the arc sine of x."); +static Py_complex asinh_special_values[7][7] = { + {{-INF,-P14},{-INF,-0.},{-INF,-0.},{-INF,0.},{-INF,0.},{-INF,P14},{-INF,N}}, + {{-INF,-P12},{U,U}, {U,U}, {U,U}, {U,U}, {-INF,P12},{N,N}}, + {{-INF,-P12},{U,U}, {-0.,-0.}, {-0.,0.}, {U,U}, {-INF,P12},{N,N}}, + {{INF,-P12}, {U,U}, {0.,-0.}, {0.,0.}, {U,U}, {INF,P12}, {N,N}}, + {{INF,-P12}, {U,U}, {U,U}, {U,U}, {U,U}, {INF,P12}, {N,N}}, + {{INF,-P14}, {INF,-0.}, {INF,-0.}, {INF,0.}, {INF,0.}, {INF,P14}, {INF,N}}, + {{INF,N}, {N,N}, {N,-0.}, {N,0.}, {N,N}, {INF,N}, {N,N}} +}; + 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; + + SPECIAL_VALUE(z, asinh_special_values); + + if (fabs(z.real) > CM_LARGE_DOUBLE || fabs(z.imag) > CM_LARGE_DOUBLE) { + if (z.imag >= 0.) { + r.real = copysign(log(hypot(z.real/2., z.imag/2.)) + + M_LN2*2., z.real); + } else { + r.real = -copysign(log(hypot(z.real/2., z.imag/2.)) + + M_LN2*2., -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); + } + errno = 0; + return r; } PyDoc_STRVAR(c_asinh_doc, @@ -86,9 +284,37 @@ PyDoc_STRVAR(c_asinh_doc, 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; +} + +/* Windows screws up atan2 for inf and nan */ +static double +c_atan2(Py_complex z) +{ + if (Py_IS_NAN(z.real) || Py_IS_NAN(z.imag)) + return Py_NAN; + if (Py_IS_INFINITY(z.imag)) { + if (Py_IS_INFINITY(z.real)) { + if (copysign(1., z.real) == 1.) + /* atan2(+-inf, +inf) == +-pi/4 */ + return copysign(0.25*Py_MATH_PI, z.imag); + else + /* atan2(+-inf, -inf) == +-pi*3/4 */ + return copysign(0.75*Py_MATH_PI, z.imag); + } + /* atan2(+-inf, x) == +-pi/2 for finite x */ + return copysign(0.5*Py_MATH_PI, z.imag); + } + return atan2(z.imag, z.real); } PyDoc_STRVAR(c_atan_doc, @@ -97,10 +323,61 @@ PyDoc_STRVAR(c_atan_doc, "Return the arc tangent of x."); +static Py_complex atanh_special_values[7][7] = { + {{-0.,-P12},{-0.,-P12},{-0.,-P12},{-0.,P12},{-0.,P12},{-0.,P12},{-0.,N}}, + {{-0.,-P12},{U,U}, {U,U}, {U,U}, {U,U}, {-0.,P12},{N,N}}, + {{-0.,-P12},{U,U}, {-0.,-0.}, {-0.,0.}, {U,U}, {-0.,P12},{-0.,N}}, + {{0.,-P12}, {U,U}, {0.,-0.}, {0.,0.}, {U,U}, {0.,P12}, {0.,N}}, + {{0.,-P12}, {U,U}, {U,U}, {U,U}, {U,U}, {0.,P12}, {N,N}}, + {{0.,-P12}, {0.,-P12}, {0.,-P12}, {0.,P12}, {0.,P12}, {0.,P12}, {0.,N}}, + {{0.,-P12}, {N,N}, {N,N}, {N,N}, {N,N}, {0.,P12}, {N,N}} +}; + 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, h; + + SPECIAL_VALUE(z, atanh_special_values); + + /* Reduce to case where z.real >= 0., using atanh(z) = -atanh(-z). */ + if (z.real < 0.) { + return c_neg(c_atanh(c_neg(z))); + } + + ay = fabs(z.imag); + if (z.real > CM_SQRT_LARGE_DOUBLE || ay > CM_SQRT_LARGE_DOUBLE) { + /* + if abs(z) is large then we use the approximation + atanh(z) ~ 1/z +/- i*pi/2 (+/- depending on the sign + of z.imag) + */ + h = hypot(z.real/2., z.imag/2.); /* safe from overflow */ + r.real = z.real/4./h/h; + /* 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(Py_MATH_PI/2., -z.imag); + errno = 0; + } else if (z.real == 1. && ay < CM_SQRT_DBL_MIN) { + /* C99 standard says: atanh(1+/-0.) should be inf +/- 0i */ + if (ay == 0.) { + r.real = INF; + r.imag = z.imag; + errno = EDOM; + } else { + r.real = -log(sqrt(ay)/sqrt(hypot(ay, 2.))); + r.imag = copysign(atan2(2., -ay)/2, z.imag); + errno = 0; + } + } 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.; + errno = 0; + } + return r; } PyDoc_STRVAR(c_atanh_doc, @@ -110,11 +387,13 @@ PyDoc_STRVAR(c_atanh_doc, 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; } @@ -124,12 +403,64 @@ PyDoc_STRVAR(c_cos_doc, "Return the cosine of x."); +/* cosh(infinity + i*y) needs to be dealt with specially */ +static Py_complex cosh_special_values[7][7] = { + {{INF,N},{U,U},{INF,0.}, {INF,-0.},{U,U},{INF,N},{INF,N}}, + {{N,N}, {U,U},{U,U}, {U,U}, {U,U},{N,N}, {N,N}}, + {{N,0.}, {U,U},{1.,0.}, {1.,-0.}, {U,U},{N,0.}, {N,0.}}, + {{N,0.}, {U,U},{1.,-0.}, {1.,0.}, {U,U},{N,0.}, {N,0.}}, + {{N,N}, {U,U},{U,U}, {U,U}, {U,U},{N,N}, {N,N}}, + {{INF,N},{U,U},{INF,-0.},{INF,0.}, {U,U},{INF,N},{INF,N}}, + {{N,N}, {N,N},{N,0.}, {N,0.}, {N,N},{N,N}, {N,N}} +}; + 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; + + /* special treatment for cosh(+/-inf + iy) if y is not a NaN */ + if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) { + if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag) && + (z.imag != 0.)) { + if (z.real > 0) { + r.real = copysign(INF, cos(z.imag)); + r.imag = copysign(INF, sin(z.imag)); + } + else { + r.real = copysign(INF, cos(z.imag)); + r.imag = -copysign(INF, sin(z.imag)); + } + } + else { + r = cosh_special_values[special_type(z.real)] + [special_type(z.imag)]; + } + /* need to set errno = EDOM if y is +/- infinity and x is not + a NaN */ + if (Py_IS_INFINITY(z.imag) && !Py_IS_NAN(z.real)) + errno = EDOM; + else + errno = 0; + return r; + } + + if (fabs(z.real) > CM_LOG_LARGE_DOUBLE) { + /* deal correctly with cases where cosh(z.real) overflows but + cosh(z) does not. */ + x_minus_one = z.real - copysign(1., z.real); + r.real = cos(z.imag) * cosh(x_minus_one) * Py_MATH_E; + r.imag = sin(z.imag) * sinh(x_minus_one) * Py_MATH_E; + } else { + r.real = cos(z.imag) * cosh(z.real); + r.imag = sin(z.imag) * sinh(z.real); + } + /* detect overflow, and set errno accordingly */ + if (Py_IS_INFINITY(r.real) || Py_IS_INFINITY(r.imag)) + errno = ERANGE; + else + errno = 0; return r; } @@ -139,13 +470,65 @@ PyDoc_STRVAR(c_cosh_doc, "Return the hyperbolic cosine of x."); +/* exp(infinity + i*y) and exp(-infinity + i*y) need special treatment for + finite y */ +static Py_complex exp_special_values[7][7] = { + {{0.,0.},{U,U},{0.,-0.}, {0.,0.}, {U,U},{0.,0.},{0.,0.}}, + {{N,N}, {U,U},{U,U}, {U,U}, {U,U},{N,N}, {N,N}}, + {{N,N}, {U,U},{1.,-0.}, {1.,0.}, {U,U},{N,N}, {N,N}}, + {{N,N}, {U,U},{1.,-0.}, {1.,0.}, {U,U},{N,N}, {N,N}}, + {{N,N}, {U,U},{U,U}, {U,U}, {U,U},{N,N}, {N,N}}, + {{INF,N},{U,U},{INF,-0.},{INF,0.},{U,U},{INF,N},{INF,N}}, + {{N,N}, {N,N},{N,-0.}, {N,0.}, {N,N},{N,N}, {N,N}} +}; + 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 (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) { + if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag) + && (z.imag != 0.)) { + if (z.real > 0) { + r.real = copysign(INF, cos(z.imag)); + r.imag = copysign(INF, sin(z.imag)); + } + else { + r.real = copysign(0., cos(z.imag)); + r.imag = copysign(0., sin(z.imag)); + } + } + else { + r = exp_special_values[special_type(z.real)] + [special_type(z.imag)]; + } + /* need to set errno = EDOM if y is +/- infinity and x is not + a NaN and not -infinity */ + if (Py_IS_INFINITY(z.imag) && + (Py_IS_FINITE(z.real) || + (Py_IS_INFINITY(z.real) && z.real > 0))) + errno = EDOM; + else + errno = 0; + return r; + } + + if (z.real > CM_LOG_LARGE_DOUBLE) { + l = exp(z.real-1.); + r.real = l*cos(z.imag)*Py_MATH_E; + r.imag = l*sin(z.imag)*Py_MATH_E; + } else { + l = exp(z.real); + r.real = l*cos(z.imag); + r.imag = l*sin(z.imag); + } + /* detect overflow, and set errno accordingly */ + if (Py_IS_INFINITY(r.real) || Py_IS_INFINITY(r.imag)) + errno = ERANGE; + else + errno = 0; return r; } @@ -155,24 +538,97 @@ PyDoc_STRVAR(c_exp_doc, "Return the exponential value e**x."); +static Py_complex log_special_values[7][7] = { + {{INF,-P34},{INF,-P}, {INF,-P}, {INF,P}, {INF,P}, {INF,P34}, {INF,N}}, + {{INF,-P12},{U,U}, {U,U}, {U,U}, {U,U}, {INF,P12}, {N,N}}, + {{INF,-P12},{U,U}, {-INF,-P}, {-INF,P}, {U,U}, {INF,P12}, {N,N}}, + {{INF,-P12},{U,U}, {-INF,-0.},{-INF,0.},{U,U}, {INF,P12}, {N,N}}, + {{INF,-P12},{U,U}, {U,U}, {U,U}, {U,U}, {INF,P12}, {N,N}}, + {{INF,-P14},{INF,-0.},{INF,-0.}, {INF,0.}, {INF,0.},{INF,P14}, {INF,N}}, + {{INF,N}, {N,N}, {N,N}, {N,N}, {N,N}, {INF,N}, {N,N}} +}; + 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 absolute error + in the log. This can be solved by rescaling z by a suitable power + of 2. + + (2) the absolute value of z is greater than DBL_MAX (e.g. when both + z.real and z.imag are within a factor of 1/sqrt(2) of DBL_MAX) + Again, rescaling solves this. + + (3) the absolute value of z is close to 1. In this case it's + difficult to achieve good accuracy, at least in part because 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. + + (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. + + */ + 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, an, h; + + SPECIAL_VALUE(z, log_special_values); + + ax = fabs(z.real); + ay = fabs(z.imag); + + if (ax > CM_LARGE_DOUBLE || ay > CM_LARGE_DOUBLE) { + r.real = log(hypot(ax/2., ay/2.)) + M_LN2; + } else if (ax < DBL_MIN && ay < DBL_MIN) { + if (ax > 0. || ay > 0.) { + /* catch cases where hypot(ax, ay) is subnormal */ + r.real = log(hypot(ldexp(ax, DBL_MANT_DIG), + ldexp(ay, DBL_MANT_DIG))) - DBL_MANT_DIG*M_LN2; + } + else { + /* log(+/-0. +/- 0i) */ + r.real = -INF; + r.imag = atan2(z.imag, z.real); + errno = EDOM; + return r; + } + } else { + h = hypot(ax, ay); + if (0.71 <= h && h <= 1.73) { + am = ax > ay ? ax : ay; /* max(ax, ay) */ + an = ax > ay ? ay : ax; /* min(ax, ay) */ + r.real = log1p((am-1)*(am+1)+an*an)/2.; + } else { + r.real = log(h); + } + } + r.imag = atan2(z.imag, z.real); + errno = 0; 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); + int errno_save; + + r = c_log(z); + errno_save = errno; /* just in case the divisions affect errno */ + r.real = r.real / M_LN10; + r.imag = r.imag / M_LN10; + errno = errno_save; return r; } @@ -182,23 +638,16 @@ PyDoc_STRVAR(c_log10_doc, "Return the base-10 logarithm of x."); -/* internal function not available from Python */ -static Py_complex -c_prodi(Py_complex x) -{ - Py_complex r; - r.real = -x.imag; - r.imag = x.real; - return r; -} - - static Py_complex -c_sin(Py_complex x) +c_sin(Py_complex z) { - Py_complex r; - r.real = sin(x.real) * cosh(x.imag); - r.imag = cos(x.real) * sinh(x.imag); + /* 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; } @@ -208,12 +657,63 @@ PyDoc_STRVAR(c_sin_doc, "Return the sine of x."); +/* sinh(infinity + i*y) needs to be dealt with specially */ +static Py_complex sinh_special_values[7][7] = { + {{INF,N},{U,U},{-INF,-0.},{-INF,0.},{U,U},{INF,N},{INF,N}}, + {{N,N}, {U,U},{U,U}, {U,U}, {U,U},{N,N}, {N,N}}, + {{0.,N}, {U,U},{-0.,-0.}, {-0.,0.}, {U,U},{0.,N}, {0.,N}}, + {{0.,N}, {U,U},{0.,-0.}, {0.,0.}, {U,U},{0.,N}, {0.,N}}, + {{N,N}, {U,U},{U,U}, {U,U}, {U,U},{N,N}, {N,N}}, + {{INF,N},{U,U},{INF,-0.}, {INF,0.}, {U,U},{INF,N},{INF,N}}, + {{N,N}, {N,N},{N,-0.}, {N,0.}, {N,N},{N,N}, {N,N}} +}; + 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; + + /* special treatment for sinh(+/-inf + iy) if y is finite and + nonzero */ + if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) { + if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag) + && (z.imag != 0.)) { + if (z.real > 0) { + r.real = copysign(INF, cos(z.imag)); + r.imag = copysign(INF, sin(z.imag)); + } + else { + r.real = -copysign(INF, cos(z.imag)); + r.imag = copysign(INF, sin(z.imag)); + } + } + else { + r = sinh_special_values[special_type(z.real)] + [special_type(z.imag)]; + } + /* need to set errno = EDOM if y is +/- infinity and x is not + a NaN */ + if (Py_IS_INFINITY(z.imag) && !Py_IS_NAN(z.real)) + errno = EDOM; + else + errno = 0; + return r; + } + + if (fabs(z.real) > CM_LOG_LARGE_DOUBLE) { + x_minus_one = z.real - copysign(1., z.real); + r.real = cos(z.imag) * sinh(x_minus_one) * Py_MATH_E; + r.imag = sin(z.imag) * cosh(x_minus_one) * Py_MATH_E; + } else { + r.real = cos(z.imag) * sinh(z.real); + r.imag = sin(z.imag) * cosh(z.real); + } + /* detect overflow, and set errno accordingly */ + if (Py_IS_INFINITY(r.real) || Py_IS_INFINITY(r.imag)) + errno = ERANGE; + else + errno = 0; return r; } @@ -223,29 +723,80 @@ PyDoc_STRVAR(c_sinh_doc, "Return the hyperbolic sine of x."); +static Py_complex sqrt_special_values[7][7] = { + {{INF,-INF},{0.,-INF},{0.,-INF},{0.,INF},{0.,INF},{INF,INF},{N,INF}}, + {{INF,-INF},{U,U}, {U,U}, {U,U}, {U,U}, {INF,INF},{N,N}}, + {{INF,-INF},{U,U}, {0.,-0.}, {0.,0.}, {U,U}, {INF,INF},{N,N}}, + {{INF,-INF},{U,U}, {0.,-0.}, {0.,0.}, {U,U}, {INF,INF},{N,N}}, + {{INF,-INF},{U,U}, {U,U}, {U,U}, {U,U}, {INF,INF},{N,N}}, + {{INF,-INF},{INF,-0.},{INF,-0.},{INF,0.},{INF,0.},{INF,INF},{INF,N}}, + {{INF,-INF},{N,N}, {N,N}, {N,N}, {N,N}, {INF,INF},{N,N}} +}; + 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/divisions, 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; + + SPECIAL_VALUE(z, sqrt_special_values); + + 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 < DBL_MIN && ay < DBL_MIN && (ax > 0. || ay > 0.)) { + /* here we catch cases where hypot(ax, ay) is subnormal */ + ax = ldexp(ax, CM_SCALE_UP); + s = ldexp(sqrt(ax + hypot(ax, ldexp(ay, CM_SCALE_UP))), + CM_SCALE_DOWN); + } 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); } + errno = 0; return r; } @@ -256,23 +807,15 @@ PyDoc_STRVAR(c_sqrt_doc, 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; } @@ -282,24 +825,78 @@ PyDoc_STRVAR(c_tan_doc, "Return the tangent of x."); +/* tanh(infinity + i*y) needs to be dealt with specially */ +static Py_complex tanh_special_values[7][7] = { + {{-1.,0.},{U,U},{-1.,-0.},{-1.,0.},{U,U},{-1.,0.},{-1.,0.}}, + {{N,N}, {U,U},{U,U}, {U,U}, {U,U},{N,N}, {N,N}}, + {{N,N}, {U,U},{-0.,-0.},{-0.,0.},{U,U},{N,N}, {N,N}}, + {{N,N}, {U,U},{0.,-0.}, {0.,0.}, {U,U},{N,N}, {N,N}}, + {{N,N}, {U,U},{U,U}, {U,U}, {U,U},{N,N}, {N,N}}, + {{1.,0.}, {U,U},{1.,-0.}, {1.,0.}, {U,U},{1.,0.}, {1.,0.}}, + {{N,N}, {N,N},{N,-0.}, {N,0.}, {N,N},{N,N}, {N,N}} +}; + 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; + + /* special treatment for tanh(+/-inf + iy) if y is finite and + nonzero */ + if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) { + if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag) + && (z.imag != 0.)) { + if (z.real > 0) { + r.real = 1.0; + r.imag = copysign(0., + 2.*sin(z.imag)*cos(z.imag)); + } + else { + r.real = -1.0; + r.imag = copysign(0., + 2.*sin(z.imag)*cos(z.imag)); + } + } + else { + r = tanh_special_values[special_type(z.real)] + [special_type(z.imag)]; + } + /* need to set errno = EDOM if z.imag is +/-infinity and + z.real is finite */ + if (Py_IS_INFINITY(z.imag) && Py_IS_FINITE(z.real)) + errno = EDOM; + else + errno = 0; + return r; + } + + /* danger of overflow in 2.*z.imag !*/ + if (fabs(z.real) > CM_LOG_LARGE_DOUBLE) { + r.real = copysign(1., z.real); + r.imag = 4.*sin(z.imag)*cos(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; + } + errno = 0; return r; } @@ -308,6 +905,7 @@ PyDoc_STRVAR(c_tanh_doc, "\n" "Return the hyperbolic tangent of x."); + static PyObject * cmath_log(PyObject *self, PyObject *args) { @@ -325,7 +923,6 @@ cmath_log(PyObject *self, PyObject *args) PyFPE_END_PROTECT(x) if (errno != 0) return math_error(); - Py_ADJUST_ERANGE2(x.real, x.imag); return PyComplex_FromCComplex(x); } @@ -351,18 +948,24 @@ math_error(void) static PyObject * math_1(PyObject *args, Py_complex (*func)(Py_complex)) { - Py_complex x; + Py_complex x,r ; if (!PyArg_ParseTuple(args, "D", &x)) return NULL; errno = 0; - PyFPE_START_PROTECT("complex function", return 0) - x = (*func)(x); - PyFPE_END_PROTECT(x) - Py_ADJUST_ERANGE2(x.real, x.imag); - if (errno != 0) - return math_error(); - else - return PyComplex_FromCComplex(x); + PyFPE_START_PROTECT("complex function", return 0); + r = (*func)(x); + PyFPE_END_PROTECT(r); + if (errno == EDOM) { + PyErr_SetString(PyExc_ValueError, "math domain error"); + return NULL; + } + else if (errno == ERANGE) { + PyErr_SetString(PyExc_OverflowError, "math range error"); + return NULL; + } + else { + return PyComplex_FromCComplex(r); + } } #define FUNC1(stubname, func) \ @@ -386,6 +989,151 @@ FUNC1(cmath_sqrt, c_sqrt) FUNC1(cmath_tan, c_tan) FUNC1(cmath_tanh, c_tanh) +static PyObject * +cmath_phase(PyObject *self, PyObject *args) +{ + Py_complex z; + double phi; + if (!PyArg_ParseTuple(args, "D:phase", &z)) + return NULL; + errno = 0; + PyFPE_START_PROTECT("arg function", return 0) + phi = c_atan2(z); + PyFPE_END_PROTECT(r) + if (errno != 0) + return math_error(); + else + return PyFloat_FromDouble(phi); +} + +PyDoc_STRVAR(cmath_phase_doc, +"phase(z) -> float\n\n\ +Return argument, also known as the phase angle, of a complex."); + +static PyObject * +cmath_polar(PyObject *self, PyObject *args) +{ + Py_complex z; + double r, phi; + if (!PyArg_ParseTuple(args, "D:polar", &z)) + return NULL; + PyFPE_START_PROTECT("polar function", return 0) + phi = c_atan2(z); /* should not cause any exception */ + r = c_abs(z); /* sets errno to ERANGE on overflow; otherwise 0 */ + PyFPE_END_PROTECT(r) + if (errno != 0) + return math_error(); + else + return Py_BuildValue("dd", r, phi); +} + +PyDoc_STRVAR(cmath_polar_doc, +"polar(z) -> r: float, phi: float\n\n\ +Convert a complex from rectangular coordinates to polar coordinates. r is\n\ +the distance from 0 and phi the phase angle."); + +/* + rect() isn't covered by the C99 standard, but it's not too hard to + figure out 'spirit of C99' rules for special value handing: + + rect(x, t) should behave like exp(log(x) + it) for positive-signed x + rect(x, t) should behave like -exp(log(-x) + it) for negative-signed x + rect(nan, t) should behave like exp(nan + it), except that rect(nan, 0) + gives nan +- i0 with the sign of the imaginary part unspecified. + +*/ + +static Py_complex rect_special_values[7][7] = { + {{INF,N},{U,U},{-INF,0.},{-INF,-0.},{U,U},{INF,N},{INF,N}}, + {{N,N}, {U,U},{U,U}, {U,U}, {U,U},{N,N}, {N,N}}, + {{0.,0.},{U,U},{-0.,0.}, {-0.,-0.}, {U,U},{0.,0.},{0.,0.}}, + {{0.,0.},{U,U},{0.,-0.}, {0.,0.}, {U,U},{0.,0.},{0.,0.}}, + {{N,N}, {U,U},{U,U}, {U,U}, {U,U},{N,N}, {N,N}}, + {{INF,N},{U,U},{INF,-0.},{INF,0.}, {U,U},{INF,N},{INF,N}}, + {{N,N}, {N,N},{N,0.}, {N,0.}, {N,N},{N,N}, {N,N}} +}; + +static PyObject * +cmath_rect(PyObject *self, PyObject *args) +{ + Py_complex z; + double r, phi; + if (!PyArg_ParseTuple(args, "dd:rect", &r, &phi)) + return NULL; + errno = 0; + PyFPE_START_PROTECT("rect function", return 0) + + /* deal with special values */ + if (!Py_IS_FINITE(r) || !Py_IS_FINITE(phi)) { + /* if r is +/-infinity and phi is finite but nonzero then + result is (+-INF +-INF i), but we need to compute cos(phi) + and sin(phi) to figure out the signs. */ + if (Py_IS_INFINITY(r) && (Py_IS_FINITE(phi) + && (phi != 0.))) { + if (r > 0) { + z.real = copysign(INF, cos(phi)); + z.imag = copysign(INF, sin(phi)); + } + else { + z.real = -copysign(INF, cos(phi)); + z.imag = -copysign(INF, sin(phi)); + } + } + else { + z = rect_special_values[special_type(r)] + [special_type(phi)]; + } + /* need to set errno = EDOM if r is a nonzero number and phi + is infinite */ + if (r != 0. && !Py_IS_NAN(r) && Py_IS_INFINITY(phi)) + errno = EDOM; + else + errno = 0; + } + else { + z.real = r * cos(phi); + z.imag = r * sin(phi); + errno = 0; + } + + PyFPE_END_PROTECT(z) + if (errno != 0) + return math_error(); + else + return PyComplex_FromCComplex(z); +} + +PyDoc_STRVAR(cmath_rect_doc, +"rect(r, phi) -> z: complex\n\n\ +Convert from polar coordinates to rectangular coordinates."); + +static PyObject * +cmath_isnan(PyObject *self, PyObject *args) +{ + Py_complex z; + if (!PyArg_ParseTuple(args, "D:isnan", &z)) + return NULL; + return PyBool_FromLong(Py_IS_NAN(z.real) || Py_IS_NAN(z.imag)); +} + +PyDoc_STRVAR(cmath_isnan_doc, +"isnan(z) -> bool\n\ +Checks if the real or imaginary part of z not a number (NaN)"); + +static PyObject * +cmath_isinf(PyObject *self, PyObject *args) +{ + Py_complex z; + if (!PyArg_ParseTuple(args, "D:isnan", &z)) + return NULL; + return PyBool_FromLong(Py_IS_INFINITY(z.real) || + Py_IS_INFINITY(z.imag)); +} + +PyDoc_STRVAR(cmath_isinf_doc, +"isinf(z) -> bool\n\ +Checks if the real or imaginary part of z is infinite."); + PyDoc_STRVAR(module_doc, "This module is always available. It provides access to mathematical\n" @@ -401,8 +1149,13 @@ static PyMethodDef cmath_methods[] = { {"cos", cmath_cos, METH_VARARGS, c_cos_doc}, {"cosh", cmath_cosh, METH_VARARGS, c_cosh_doc}, {"exp", cmath_exp, METH_VARARGS, c_exp_doc}, + {"isinf", cmath_isinf, METH_VARARGS, cmath_isinf_doc}, + {"isnan", cmath_isnan, METH_VARARGS, cmath_isnan_doc}, {"log", cmath_log, METH_VARARGS, cmath_log_doc}, {"log10", cmath_log10, METH_VARARGS, c_log10_doc}, + {"phase", cmath_phase, METH_VARARGS, cmath_phase_doc}, + {"polar", cmath_polar, METH_VARARGS, cmath_polar_doc}, + {"rect", cmath_rect, METH_VARARGS, cmath_rect_doc}, {"sin", cmath_sin, METH_VARARGS, c_sin_doc}, {"sinh", cmath_sinh, METH_VARARGS, c_sinh_doc}, {"sqrt", cmath_sqrt, METH_VARARGS, c_sqrt_doc}, @@ -421,6 +1174,6 @@ initcmath(void) return; PyModule_AddObject(m, "pi", - PyFloat_FromDouble(atan(1.0) * 4.0)); - PyModule_AddObject(m, "e", PyFloat_FromDouble(exp(1.0))); + PyFloat_FromDouble(Py_MATH_PI)); + PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E)); } diff --git a/Modules/mathmodule.c b/Modules/mathmodule.c index cf2bf64..8c48316 100644 --- a/Modules/mathmodule.c +++ b/Modules/mathmodule.c @@ -1,17 +1,60 @@ /* Math module -- standard C math library functions, pi and e */ +/* Here are some comments from Tim Peters, extracted from the + discussion attached to http://bugs.python.org/issue1640. They + describe the general aims of the math module with respect to + special values, IEEE-754 floating-point exceptions, and Python + exceptions. + +These are the "spirit of 754" rules: + +1. If the mathematical result is a real number, but of magnitude too +large to approximate by a machine float, overflow is signaled and the +result is an infinity (with the appropriate sign). + +2. If the mathematical result is a real number, but of magnitude too +small to approximate by a machine float, underflow is signaled and the +result is a zero (with the appropriate sign). + +3. At a singularity (a value x such that the limit of f(y) as y +approaches x exists and is an infinity), "divide by zero" is signaled +and the result is an infinity (with the appropriate sign). This is +complicated a little by that the left-side and right-side limits may +not be the same; e.g., 1/x approaches +inf or -inf as x approaches 0 +from the positive or negative directions. In that specific case, the +sign of the zero determines the result of 1/0. + +4. At a point where a function has no defined result in the extended +reals (i.e., the reals plus an infinity or two), invalid operation is +signaled and a NaN is returned. + +And these are what Python has historically /tried/ to do (but not +always successfully, as platform libm behavior varies a lot): + +For #1, raise OverflowError. + +For #2, return a zero (with the appropriate sign if that happens by +accident ;-)). + +For #3 and #4, raise ValueError. It may have made sense to raise +Python's ZeroDivisionError in #3, but historically that's only been +raised for division by zero and mod by zero. + +*/ + +/* + In general, on an IEEE-754 platform the aim is to follow the C99 + standard, including Annex 'F', whenever possible. Where the + standard recommends raising the 'divide-by-zero' or 'invalid' + floating-point exceptions, Python should raise a ValueError. Where + the standard recommends raising 'overflow', Python should raise an + OverflowError. In all other circumstances a value should be + returned. + */ + #include "Python.h" #include "longintrepr.h" /* just for SHIFT */ -#ifndef _MSC_VER -#ifndef __STDC__ -extern double fmod (double, double); -extern double frexp (double, int *); -extern double ldexp (double, int); -extern double modf (double, double *); -#endif /* __STDC__ */ -#endif /* _MSC_VER */ - #ifdef _OSF_SOURCE /* OSF1 5.1 doesn't make this available with XOPEN_SOURCE_EXTENDED defined */ extern double copysign(double, double); @@ -52,41 +95,111 @@ is_error(double x) return result; } +/* + math_1 is used to wrap a libm function f that takes a double + arguments and returns a double. + + The error reporting follows these rules, which are designed to do + the right thing on C89/C99 platforms and IEEE 754/non IEEE 754 + platforms. + + - a NaN result from non-NaN inputs causes ValueError to be raised + - an infinite result from finite inputs causes OverflowError to be + raised if can_overflow is 1, or raises ValueError if can_overflow + is 0. + - if the result is finite and errno == EDOM then ValueError is + raised + - if the result is finite and nonzero and errno == ERANGE then + OverflowError is raised + + The last rule is used to catch overflow on platforms which follow + C89 but for which HUGE_VAL is not an infinity. + + For the majority of one-argument functions these rules are enough + to ensure that Python's functions behave as specified in 'Annex F' + of the C99 standard, with the 'invalid' and 'divide-by-zero' + floating-point exceptions mapping to Python's ValueError and the + 'overflow' floating-point exception mapping to OverflowError. + math_1 only works for functions that don't have singularities *and* + the possibility of overflow; fortunately, that covers everything we + care about right now. +*/ + static PyObject * math_1_to_whatever(PyObject *arg, double (*func) (double), - PyObject *(*from_double_func) (double)) + PyObject *(*from_double_func) (double), + int can_overflow) { - double x = PyFloat_AsDouble(arg); + double x, r; + x = PyFloat_AsDouble(arg); if (x == -1.0 && PyErr_Occurred()) return NULL; errno = 0; - PyFPE_START_PROTECT("in math_1", return 0) - x = (*func)(x); - PyFPE_END_PROTECT(x) - Py_SET_ERRNO_ON_MATH_ERROR(x); - if (errno && is_error(x)) + PyFPE_START_PROTECT("in math_1", return 0); + r = (*func)(x); + PyFPE_END_PROTECT(r); + if (Py_IS_NAN(r)) { + if (!Py_IS_NAN(x)) + errno = EDOM; + else + errno = 0; + } + else if (Py_IS_INFINITY(r)) { + if (Py_IS_FINITE(x)) + errno = can_overflow ? ERANGE : EDOM; + else + errno = 0; + } + if (errno && is_error(r)) return NULL; else - return (*from_double_func)(x); + return (*from_double_func)(r); } +/* + math_2 is used to wrap a libm function f that takes two double + arguments and returns a double. + + The error reporting follows these rules, which are designed to do + the right thing on C89/C99 platforms and IEEE 754/non IEEE 754 + platforms. + + - a NaN result from non-NaN inputs causes ValueError to be raised + - an infinite result from finite inputs causes OverflowError to be + raised. + - if the result is finite and errno == EDOM then ValueError is + raised + - if the result is finite and nonzero and errno == ERANGE then + OverflowError is raised + + The last rule is used to catch overflow on platforms which follow + C89 but for which HUGE_VAL is not an infinity. + + For most two-argument functions (copysign, fmod, hypot, atan2) + these rules are enough to ensure that Python's functions behave as + specified in 'Annex F' of the C99 standard, with the 'invalid' and + 'divide-by-zero' floating-point exceptions mapping to Python's + ValueError and the 'overflow' floating-point exception mapping to + OverflowError. +*/ + static PyObject * -math_1(PyObject *arg, double (*func) (double)) +math_1(PyObject *arg, double (*func) (double), int can_overflow) { - return math_1_to_whatever(arg, func, PyFloat_FromDouble); + return math_1_to_whatever(arg, func, PyFloat_FromDouble, can_overflow); } static PyObject * -math_1_to_int(PyObject *arg, double (*func) (double)) +math_1_to_int(PyObject *arg, double (*func) (double), int can_overflow) { - return math_1_to_whatever(arg, func, PyLong_FromDouble); + return math_1_to_whatever(arg, func, PyLong_FromDouble, can_overflow); } static PyObject * math_2(PyObject *args, double (*func) (double, double), char *funcname) { PyObject *ox, *oy; - double x, y; + double x, y, r; if (! PyArg_UnpackTuple(args, funcname, 2, 2, &ox, &oy)) return NULL; x = PyFloat_AsDouble(ox); @@ -94,19 +207,30 @@ math_2(PyObject *args, double (*func) (double, double), char *funcname) if ((x == -1.0 || y == -1.0) && PyErr_Occurred()) return NULL; errno = 0; - PyFPE_START_PROTECT("in math_2", return 0) - x = (*func)(x, y); - PyFPE_END_PROTECT(x) - Py_SET_ERRNO_ON_MATH_ERROR(x); - if (errno && is_error(x)) + PyFPE_START_PROTECT("in math_2", return 0); + r = (*func)(x, y); + PyFPE_END_PROTECT(r); + if (Py_IS_NAN(r)) { + if (!Py_IS_NAN(x) && !Py_IS_NAN(y)) + errno = EDOM; + else + errno = 0; + } + else if (Py_IS_INFINITY(r)) { + if (Py_IS_FINITE(x) && Py_IS_FINITE(y)) + errno = ERANGE; + else + errno = 0; + } + if (errno && is_error(r)) return NULL; else - return PyFloat_FromDouble(x); + return PyFloat_FromDouble(r); } -#define FUNC1(funcname, func, docstring) \ +#define FUNC1(funcname, func, can_overflow, docstring) \ static PyObject * math_##funcname(PyObject *self, PyObject *args) { \ - return math_1(args, func); \ + return math_1(args, func, can_overflow); \ }\ PyDoc_STRVAR(math_##funcname##_doc, docstring); @@ -116,15 +240,21 @@ math_2(PyObject *args, double (*func) (double, double), char *funcname) }\ PyDoc_STRVAR(math_##funcname##_doc, docstring); -FUNC1(acos, acos, +FUNC1(acos, acos, 0, "acos(x)\n\nReturn the arc cosine (measured in radians) of x.") -FUNC1(asin, asin, +FUNC1(acosh, acosh, 0, + "acosh(x)\n\nReturn the hyperbolic arc cosine (measured in radians) of x.") +FUNC1(asin, asin, 0, "asin(x)\n\nReturn the arc sine (measured in radians) of x.") -FUNC1(atan, atan, +FUNC1(asinh, asinh, 0, + "asinh(x)\n\nReturn the hyperbolic arc sine (measured in radians) of x.") +FUNC1(atan, atan, 0, "atan(x)\n\nReturn the arc tangent (measured in radians) of x.") FUNC2(atan2, atan2, "atan2(y, x)\n\nReturn the arc tangent (measured in radians) of y/x.\n" "Unlike atan(y/x), the signs of both x and y are considered.") +FUNC1(atanh, atanh, 0, + "atanh(x)\n\nReturn the hyperbolic arc tangent (measured in radians) of x.") static PyObject * math_ceil(PyObject *self, PyObject *number) { static PyObject *ceil_str = NULL; @@ -138,7 +268,7 @@ static PyObject * math_ceil(PyObject *self, PyObject *number) { method = _PyType_Lookup(Py_TYPE(number), ceil_str); if (method == NULL) - return math_1_to_int(number, ceil); + return math_1_to_int(number, ceil, 0); else return PyObject_CallFunction(method, "O", number); } @@ -147,23 +277,15 @@ PyDoc_STRVAR(math_ceil_doc, "ceil(x)\n\nReturn the ceiling of x as an int.\n" "This is the smallest integral value >= x."); -FUNC1(cos, cos, +FUNC2(copysign, copysign, + "copysign(x,y)\n\nReturn x with the sign of y.") +FUNC1(cos, cos, 0, "cos(x)\n\nReturn the cosine of x (measured in radians).") -FUNC1(cosh, cosh, +FUNC1(cosh, cosh, 1, "cosh(x)\n\nReturn the hyperbolic cosine of x.") - -#ifdef MS_WINDOWS -# define copysign _copysign -# define HAVE_COPYSIGN 1 -#endif -#ifdef HAVE_COPYSIGN -FUNC2(copysign, copysign, - "copysign(x,y)\n\nReturn x with the sign of y."); -#endif - -FUNC1(exp, exp, +FUNC1(exp, exp, 1, "exp(x)\n\nReturn e raised to the power of x.") -FUNC1(fabs, fabs, +FUNC1(fabs, fabs, 0, "fabs(x)\n\nReturn the absolute value of the float x.") static PyObject * math_floor(PyObject *self, PyObject *number) { @@ -178,7 +300,7 @@ static PyObject * math_floor(PyObject *self, PyObject *number) { method = _PyType_Lookup(Py_TYPE(number), floor_str); if (method == NULL) - return math_1_to_int(number, floor); + return math_1_to_int(number, floor, 0); else return PyObject_CallFunction(method, "O", number); } @@ -187,22 +309,18 @@ PyDoc_STRVAR(math_floor_doc, "floor(x)\n\nReturn the floor of x as an int.\n" "This is the largest integral value <= x."); -FUNC2(fmod, fmod, - "fmod(x,y)\n\nReturn fmod(x, y), according to platform C." - " x % y may differ.") -FUNC2(hypot, hypot, - "hypot(x,y)\n\nReturn the Euclidean distance, sqrt(x*x + y*y).") -FUNC2(pow, pow, - "pow(x,y)\n\nReturn x**y (x to the power of y).") -FUNC1(sin, sin, +FUNC1(log1p, log1p, 1, + "log1p(x)\n\nReturn the natural logarithm of 1+x (base e).\n\ + The result is computed in a way which is accurate for x near zero.") +FUNC1(sin, sin, 0, "sin(x)\n\nReturn the sine of x (measured in radians).") -FUNC1(sinh, sinh, +FUNC1(sinh, sinh, 1, "sinh(x)\n\nReturn the hyperbolic sine of x.") -FUNC1(sqrt, sqrt, +FUNC1(sqrt, sqrt, 0, "sqrt(x)\n\nReturn the square root of x.") -FUNC1(tan, tan, +FUNC1(tan, tan, 0, "tan(x)\n\nReturn the tangent of x (measured in radians).") -FUNC1(tanh, tanh, +FUNC1(tanh, tanh, 0, "tanh(x)\n\nReturn the hyperbolic tangent of x.") static PyObject * @@ -244,13 +362,17 @@ math_frexp(PyObject *self, PyObject *arg) double x = PyFloat_AsDouble(arg); if (x == -1.0 && PyErr_Occurred()) return NULL; - errno = 0; - x = frexp(x, &i); - Py_SET_ERRNO_ON_MATH_ERROR(x); - if (errno && is_error(x)) - return NULL; - else - return Py_BuildValue("(di)", x, i); + /* deal with special cases directly, to sidestep platform + differences */ + if (Py_IS_NAN(x) || Py_IS_INFINITY(x) || !x) { + i = 0; + } + else { + PyFPE_START_PROTECT("in math_frexp", return 0); + x = frexp(x, &i); + PyFPE_END_PROTECT(x); + } + return Py_BuildValue("(di)", x, i); } PyDoc_STRVAR(math_frexp_doc, @@ -263,19 +385,24 @@ PyDoc_STRVAR(math_frexp_doc, static PyObject * math_ldexp(PyObject *self, PyObject *args) { - double x; + double x, r; int exp; if (! PyArg_ParseTuple(args, "di:ldexp", &x, &exp)) return NULL; errno = 0; - PyFPE_START_PROTECT("ldexp", return 0) - x = ldexp(x, exp); - PyFPE_END_PROTECT(x) - Py_SET_ERRNO_ON_MATH_ERROR(x); - if (errno && is_error(x)) + PyFPE_START_PROTECT("in math_ldexp", return 0) + r = ldexp(x, exp); + PyFPE_END_PROTECT(r) + if (Py_IS_FINITE(x) && Py_IS_INFINITY(r)) + errno = ERANGE; + /* Windows MSVC8 sets errno = EDOM on ldexp(NaN, i); + we unset it to avoid raising a ValueError here. */ + if (errno == EDOM) + errno = 0; + if (errno && is_error(r)) return NULL; else - return PyFloat_FromDouble(x); + return PyFloat_FromDouble(r); } PyDoc_STRVAR(math_ldexp_doc, @@ -288,12 +415,10 @@ math_modf(PyObject *self, PyObject *arg) if (x == -1.0 && PyErr_Occurred()) return NULL; errno = 0; + PyFPE_START_PROTECT("in math_modf", return 0); x = modf(x, &y); - Py_SET_ERRNO_ON_MATH_ERROR(x); - if (errno && is_error(x)) - return NULL; - else - return Py_BuildValue("(dd)", x, y); + PyFPE_END_PROTECT(x); + return Py_BuildValue("(dd)", x, y); } PyDoc_STRVAR(math_modf_doc, @@ -332,7 +457,7 @@ loghelper(PyObject* arg, double (*func)(double), char *funcname) } /* Else let libm handle it by itself. */ - return math_1(arg, func); + return math_1(arg, func, 0); } static PyObject * @@ -375,6 +500,141 @@ math_log10(PyObject *self, PyObject *arg) PyDoc_STRVAR(math_log10_doc, "log10(x) -> the base 10 logarithm of x."); +static PyObject * +math_fmod(PyObject *self, PyObject *args) +{ + PyObject *ox, *oy; + double r, x, y; + if (! PyArg_UnpackTuple(args, "fmod", 2, 2, &ox, &oy)) + return NULL; + x = PyFloat_AsDouble(ox); + y = PyFloat_AsDouble(oy); + if ((x == -1.0 || y == -1.0) && PyErr_Occurred()) + return NULL; + /* fmod(x, +/-Inf) returns x for finite x. */ + if (Py_IS_INFINITY(y) && Py_IS_FINITE(x)) + return PyFloat_FromDouble(x); + errno = 0; + PyFPE_START_PROTECT("in math_fmod", return 0); + r = fmod(x, y); + PyFPE_END_PROTECT(r); + if (Py_IS_NAN(r)) { + if (!Py_IS_NAN(x) && !Py_IS_NAN(y)) + errno = EDOM; + else + errno = 0; + } + if (errno && is_error(r)) + return NULL; + else + return PyFloat_FromDouble(r); +} + +PyDoc_STRVAR(math_fmod_doc, +"fmod(x,y)\n\nReturn fmod(x, y), according to platform C." +" x % y may differ."); + +static PyObject * +math_hypot(PyObject *self, PyObject *args) +{ + PyObject *ox, *oy; + double r, x, y; + if (! PyArg_UnpackTuple(args, "hypot", 2, 2, &ox, &oy)) + return NULL; + x = PyFloat_AsDouble(ox); + y = PyFloat_AsDouble(oy); + if ((x == -1.0 || y == -1.0) && PyErr_Occurred()) + return NULL; + /* hypot(x, +/-Inf) returns Inf, even if x is a NaN. */ + if (Py_IS_INFINITY(x)) + return PyFloat_FromDouble(fabs(x)); + if (Py_IS_INFINITY(y)) + return PyFloat_FromDouble(fabs(y)); + errno = 0; + PyFPE_START_PROTECT("in math_hypot", return 0); + r = hypot(x, y); + PyFPE_END_PROTECT(r); + if (Py_IS_NAN(r)) { + if (!Py_IS_NAN(x) && !Py_IS_NAN(y)) + errno = EDOM; + else + errno = 0; + } + else if (Py_IS_INFINITY(r)) { + if (Py_IS_FINITE(x) && Py_IS_FINITE(y)) + errno = ERANGE; + else + errno = 0; + } + if (errno && is_error(r)) + return NULL; + else + return PyFloat_FromDouble(r); +} + +PyDoc_STRVAR(math_hypot_doc, +"hypot(x,y)\n\nReturn the Euclidean distance, sqrt(x*x + y*y)."); + +/* pow can't use math_2, but needs its own wrapper: the problem is + that an infinite result can arise either as a result of overflow + (in which case OverflowError should be raised) or as a result of + e.g. 0.**-5. (for which ValueError needs to be raised.) +*/ + +static PyObject * +math_pow(PyObject *self, PyObject *args) +{ + PyObject *ox, *oy; + double r, x, y; + + if (! PyArg_UnpackTuple(args, "pow", 2, 2, &ox, &oy)) + return NULL; + x = PyFloat_AsDouble(ox); + y = PyFloat_AsDouble(oy); + if ((x == -1.0 || y == -1.0) && PyErr_Occurred()) + return NULL; + /* 1**x and x**0 return 1., even if x is a NaN or infinity. */ + if (x == 1.0 || y == 0.0) + return PyFloat_FromDouble(1.); + errno = 0; + PyFPE_START_PROTECT("in math_pow", return 0); + r = pow(x, y); + PyFPE_END_PROTECT(r); + if (Py_IS_NAN(r)) { + if (!Py_IS_NAN(x) && !Py_IS_NAN(y)) + errno = EDOM; + else + errno = 0; + } + /* an infinite result arises either from: + + (A) (+/-0.)**negative, + (B) overflow of x**y with both x and y finite (and x nonzero) + (C) (+/-inf)**positive, or + (D) x**inf with |x| > 1, or x**-inf with |x| < 1. + + In case (A) we want ValueError to be raised. In case (B) + OverflowError should be raised. In cases (C) and (D) the infinite + result should be returned. + */ + else if (Py_IS_INFINITY(r)) { + if (x == 0.) + errno = EDOM; + else if (Py_IS_FINITE(x) && Py_IS_FINITE(y)) + errno = ERANGE; + else + errno = 0; + } + + if (errno && is_error(r)) + return NULL; + else + return PyFloat_FromDouble(r); +} + +PyDoc_STRVAR(math_pow_doc, +"pow(x,y)\n\nReturn x**y (x to the power of y)."); + static const double degToRad = Py_MATH_PI / 180.0; static const double radToDeg = 180.0 / Py_MATH_PI; @@ -428,16 +688,16 @@ PyDoc_STRVAR(math_isinf_doc, "isinf(x) -> bool\n\ Checks if float x is infinite (positive or negative)"); - static PyMethodDef math_methods[] = { {"acos", math_acos, METH_O, math_acos_doc}, + {"acosh", math_acosh, METH_O, math_acosh_doc}, {"asin", math_asin, METH_O, math_asin_doc}, + {"asinh", math_asinh, METH_O, math_asinh_doc}, {"atan", math_atan, METH_O, math_atan_doc}, {"atan2", math_atan2, METH_VARARGS, math_atan2_doc}, + {"atanh", math_atanh, METH_O, math_atanh_doc}, {"ceil", math_ceil, METH_O, math_ceil_doc}, -#ifdef HAVE_COPYSIGN {"copysign", math_copysign, METH_VARARGS, math_copysign_doc}, -#endif {"cos", math_cos, METH_O, math_cos_doc}, {"cosh", math_cosh, METH_O, math_cosh_doc}, {"degrees", math_degrees, METH_O, math_degrees_doc}, @@ -451,6 +711,7 @@ static PyMethodDef math_methods[] = { {"isnan", math_isnan, METH_O, math_isnan_doc}, {"ldexp", math_ldexp, METH_VARARGS, math_ldexp_doc}, {"log", math_log, METH_VARARGS, math_log_doc}, + {"log1p", math_log1p, METH_O, math_log1p_doc}, {"log10", math_log10, METH_O, math_log10_doc}, {"modf", math_modf, METH_O, math_modf_doc}, {"pow", math_pow, METH_VARARGS, math_pow_doc}, @@ -472,27 +733,15 @@ PyDoc_STRVAR(module_doc, PyMODINIT_FUNC initmath(void) { - PyObject *m, *d, *v; + PyObject *m; m = Py_InitModule3("math", math_methods, module_doc); if (m == NULL) goto finally; - d = PyModule_GetDict(m); - if (d == NULL) - goto finally; - - if (!(v = PyFloat_FromDouble(Py_MATH_PI))) - goto finally; - if (PyDict_SetItemString(d, "pi", v) < 0) - goto finally; - Py_DECREF(v); - if (!(v = PyFloat_FromDouble(Py_MATH_E))) - goto finally; - if (PyDict_SetItemString(d, "e", v) < 0) - goto finally; - Py_DECREF(v); + PyModule_AddObject(m, "pi", PyFloat_FromDouble(Py_MATH_PI)); + PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E)); - finally: + finally: return; } |