Eigenschaftsänderungen: . ___________________________________________________________________ Name: svnmerge-integrated - /python/branches/libffi3-branch:1-61232 /python/branches/trunk-math:1-60195 + /python/branches/libffi3-branch:1-61232 /python/branches/trunk-math:1-61829 Index: Python/hypot.c =================================================================== --- Python/hypot.c (Revision 61828) +++ Python/hypot.c (Arbeitskopie) @@ -1,25 +0,0 @@ -/* hypot() replacement */ - -#include "Python.h" - -#ifndef HAVE_HYPOT -double hypot(double x, double y) -{ - double yx; - - x = fabs(x); - y = fabs(y); - if (x < y) { - double temp = x; - x = y; - y = temp; - } - if (x == 0.) - return 0.; - else { - yx = y/x; - return x*sqrt(1.+yx*yx); - } -} -#endif /* HAVE_HYPOT */ - Index: configure =================================================================== --- configure (Revision 61828) +++ configure (Arbeitskopie) @@ -1,5 +1,5 @@ #! /bin/sh -# From configure.in Revision: 61436 . +# From configure.in Revision: 61722 . # Guess values for system-dependent variables and create Makefiles. # Generated by GNU Autoconf 2.61 for python 2.6. # Index: Include/pyport.h =================================================================== --- Include/pyport.h (Revision 61828) +++ Include/pyport.h (Arbeitskopie) @@ -353,123 +353,6 @@ #define Py_SAFE_DOWNCAST(VALUE, WIDE, NARROW) (NARROW)(VALUE) #endif -/* High precision defintion of pi and e (Euler) - * The values are taken from libc6's math.h. - */ -#ifndef Py_MATH_PIl -#define Py_MATH_PIl 3.1415926535897932384626433832795029L -#endif -#ifndef Py_MATH_PI -#define Py_MATH_PI 3.14159265358979323846 -#endif - -#ifndef Py_MATH_El -#define Py_MATH_El 2.7182818284590452353602874713526625L -#endif - -#ifndef Py_MATH_E -#define Py_MATH_E 2.7182818284590452354 -#endif - -/* Py_IS_NAN(X) - * Return 1 if float or double arg is a NaN, else 0. - * Caution: - * X is evaluated more than once. - * This may not work on all platforms. Each platform has *some* - * way to spell this, though -- override in pyconfig.h if you have - * a platform where it doesn't work. - */ -#ifndef Py_IS_NAN -#ifdef HAVE_ISNAN -#define Py_IS_NAN(X) isnan(X) -#else -#define Py_IS_NAN(X) ((X) != (X)) -#endif -#endif - -/* Py_IS_INFINITY(X) - * Return 1 if float or double arg is an infinity, else 0. - * Caution: - * X is evaluated more than once. - * This implementation may set the underflow flag if |X| is very small; - * it really can't be implemented correctly (& easily) before C99. - * Override in pyconfig.h if you have a better spelling on your platform. - */ -#ifndef Py_IS_INFINITY -#ifdef HAVE_ISINF -#define Py_IS_INFINITY(X) isinf(X) -#else -#define Py_IS_INFINITY(X) ((X) && (X)*0.5 == (X)) -#endif -#endif - -/* Py_IS_FINITE(X) - * Return 1 if float or double arg is neither infinite nor NAN, else 0. - * Some compilers (e.g. VisualStudio) have intrisics for this, so a special - * macro for this particular test is useful - */ -#ifndef Py_IS_FINITE -#ifdef HAVE_FINITE -#define Py_IS_FINITE(X) finite(X) -#else -#define Py_IS_FINITE(X) (!Py_IS_INFINITY(X) && !Py_IS_NAN(X)) -#endif -#endif - -/* HUGE_VAL is supposed to expand to a positive double infinity. Python - * uses Py_HUGE_VAL instead because some platforms are broken in this - * respect. We used to embed code in pyport.h to try to worm around that, - * but different platforms are broken in conflicting ways. If you're on - * a platform where HUGE_VAL is defined incorrectly, fiddle your Python - * config to #define Py_HUGE_VAL to something that works on your platform. - */ -#ifndef Py_HUGE_VAL -#define Py_HUGE_VAL HUGE_VAL -#endif - -/* Py_NAN - * A value that evaluates to a NaN. On IEEE 754 platforms INF*0 or - * INF/INF works. Define Py_NO_NAN in pyconfig.h if your platform - * doesn't support NaNs. - */ -#if !defined(Py_NAN) && !defined(Py_NO_NAN) -#define Py_NAN (Py_HUGE_VAL * 0.) -#endif - -/* Py_OVERFLOWED(X) - * Return 1 iff a libm function overflowed. Set errno to 0 before calling - * a libm function, and invoke this macro after, passing the function - * result. - * Caution: - * This isn't reliable. C99 no longer requires libm to set errno under - * any exceptional condition, but does require +- HUGE_VAL return - * values on overflow. A 754 box *probably* maps HUGE_VAL to a - * double infinity, and we're cool if that's so, unless the input - * was an infinity and an infinity is the expected result. A C89 - * system sets errno to ERANGE, so we check for that too. We're - * out of luck if a C99 754 box doesn't map HUGE_VAL to +Inf, or - * if the returned result is a NaN, or if a C89 box returns HUGE_VAL - * in non-overflow cases. - * X is evaluated more than once. - * Some platforms have better way to spell this, so expect some #ifdef'ery. - * - * OpenBSD uses 'isinf()' because a compiler bug on that platform causes - * the longer macro version to be mis-compiled. This isn't optimal, and - * should be removed once a newer compiler is available on that platform. - * The system that had the failure was running OpenBSD 3.2 on Intel, with - * gcc 2.95.3. - * - * According to Tim's checkin, the FreeBSD systems use isinf() to work - * around a FPE bug on that platform. - */ -#if defined(__FreeBSD__) || defined(__OpenBSD__) -#define Py_OVERFLOWED(X) isinf(X) -#else -#define Py_OVERFLOWED(X) ((X) != 0.0 && (errno == ERANGE || \ - (X) == Py_HUGE_VAL || \ - (X) == -Py_HUGE_VAL)) -#endif - /* Py_SET_ERRNO_ON_MATH_ERROR(x) * If a libm function did not set errno, but it looks like the result * overflowed or not-a-number, set errno to ERANGE or EDOM. Set errno @@ -601,15 +484,6 @@ #endif /* 0 */ -/************************ - * WRAPPER FOR * - ************************/ - -#ifndef HAVE_HYPOT -extern double hypot(double, double); -#endif - - /* On 4.4BSD-descendants, ctype functions serves the whole range of * wchar_t character set rather than single byte code points only. * This characteristic can break some operations of string object Index: Include/complexobject.h =================================================================== --- Include/complexobject.h (Revision 61828) +++ Include/complexobject.h (Arbeitskopie) @@ -19,6 +19,7 @@ #define c_prod _Py_c_prod #define c_quot _Py_c_quot #define c_pow _Py_c_pow +#define c_abs _Py_c_abs PyAPI_FUNC(Py_complex) c_sum(Py_complex, Py_complex); PyAPI_FUNC(Py_complex) c_diff(Py_complex, Py_complex); @@ -26,6 +27,7 @@ PyAPI_FUNC(Py_complex) c_prod(Py_complex, Py_complex); PyAPI_FUNC(Py_complex) c_quot(Py_complex, Py_complex); PyAPI_FUNC(Py_complex) c_pow(Py_complex, Py_complex); +PyAPI_FUNC(double) c_abs(Py_complex); /* Complex object interface */ Index: Include/Python.h =================================================================== --- Include/Python.h (Revision 61828) +++ Include/Python.h (Arbeitskopie) @@ -73,6 +73,7 @@ #if defined(PYMALLOC_DEBUG) && !defined(WITH_PYMALLOC) #error "PYMALLOC_DEBUG requires WITH_PYMALLOC" #endif +#include "pymath.h" #include "pymem.h" #include "object.h" Index: Include/floatobject.h =================================================================== --- Include/floatobject.h (Revision 61828) +++ Include/floatobject.h (Arbeitskopie) @@ -21,6 +21,17 @@ #define PyFloat_Check(op) PyObject_TypeCheck(op, &PyFloat_Type) #define PyFloat_CheckExact(op) (Py_TYPE(op) == &PyFloat_Type) +#ifdef Py_NAN +#define Py_RETURN_NAN return Py_INCREF(PyFloat_NAN), PyFloat_NAN +#endif + +#define Py_RETURN_INF(sign) do \ + if (copysign(1., sign) == 1.) { \ + return Py_INCREF(PyFloat_PINF), PyFloat_PINF; \ + } else { \ + return Py_INCREF(PyFloat_NINF), PyFloat_NINF; \ + } while(0) + PyAPI_FUNC(double) PyFloat_GetMax(void); PyAPI_FUNC(double) PyFloat_GetMin(void); PyAPI_FUNC(PyObject *) PyFloat_GetInfo(void); Index: Objects/complexobject.c =================================================================== --- Objects/complexobject.c (Revision 61828) +++ Objects/complexobject.c (Arbeitskopie) @@ -183,6 +183,38 @@ } +double +c_abs(Py_complex z) +{ + /* sets errno = ERANGE on overflow; otherwise errno = 0 */ + double result; + + if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) { + /* C99 rules: if either the real or the imaginary part is an + infinity, return infinity, even if the other part is a + NaN. */ + if (Py_IS_INFINITY(z.real)) { + result = fabs(z.real); + errno = 0; + return result; + } + if (Py_IS_INFINITY(z.imag)) { + result = fabs(z.imag); + errno = 0; + return result; + } + /* either the real or imaginary part is a NaN, + and neither is infinite. Result should be NaN. */ + return Py_NAN; + } + result = hypot(z.real, z.imag); + if (!Py_IS_FINITE(result)) + errno = ERANGE; + else + errno = 0; + return result; +} + static PyObject * complex_subtype_from_c_complex(PyTypeObject *type, Py_complex cval) { @@ -325,8 +357,7 @@ if (!Py_IS_FINITE(v->cval.imag)) { if (Py_IS_NAN(v->cval.imag)) strncpy(buf, "nan*j", 6); - /* else if (copysign(1, v->cval.imag) == 1) */ - else if (v->cval.imag > 0) + else if (copysign(1, v->cval.imag) == 1) strncpy(buf, "inf*j", 6); else strncpy(buf, "-inf*j", 7); @@ -488,6 +519,7 @@ complex_div(PyComplexObject *v, PyComplexObject *w) { Py_complex quot; + PyFPE_START_PROTECT("complex_div", return 0) errno = 0; quot = c_quot(v->cval,w->cval); @@ -651,9 +683,16 @@ complex_abs(PyComplexObject *v) { double result; + PyFPE_START_PROTECT("complex_abs", return 0) - result = hypot(v->cval.real,v->cval.imag); + result = c_abs(v->cval); PyFPE_END_PROTECT(result) + + if (errno == ERANGE) { + PyErr_SetString(PyExc_OverflowError, + "absolute value too large"); + return NULL; + } return PyFloat_FromDouble(result); } @@ -782,9 +821,25 @@ return Py_BuildValue("(D)", &v->cval); } +static PyObject * +complex_is_finite(PyObject *self) +{ + Py_complex c; + c = ((PyComplexObject *)self)->cval; + return PyBool_FromLong((long)(Py_IS_FINITE(c.real) && + Py_IS_FINITE(c.imag))); +} + +PyDoc_STRVAR(complex_is_finite_doc, +"complex.is_finite() -> bool\n" +"\n" +"Returns True if the real and the imaginary part is finite."); + static PyMethodDef complex_methods[] = { {"conjugate", (PyCFunction)complex_conjugate, METH_NOARGS, complex_conjugate_doc}, + {"is_finite", (PyCFunction)complex_is_finite, METH_NOARGS, + complex_is_finite_doc}, {"__getnewargs__", (PyCFunction)complex_getnewargs, METH_NOARGS}, {NULL, NULL} /* sentinel */ }; Index: Objects/doubledigits.c =================================================================== --- Objects/doubledigits.c (Revision 61828) +++ Objects/doubledigits.c (Arbeitskopie) @@ -1,601 +0,0 @@ -/* Free-format floating point printer - * - * Based on "Floating-Point Printer Sample Code", by Robert G. Burger, - * http://www.cs.indiana.edu/~burger/fp/index.html - */ - -#include "Python.h" - -#if defined(__alpha) || defined(__i386) || defined(_M_IX86) || defined(_M_X64) || defined(_M_IA64) -#define LITTLE_ENDIAN_IEEE_DOUBLE -#elif !(defined(__ppc__) || defined(sparc) || defined(__sgi) || defined(_IBMR2) || defined(hpux)) -#error unknown machine type -#endif - -#if defined(_M_IX86) -#define UNSIGNED64 unsigned __int64 -#elif defined(__alpha) -#define UNSIGNED64 unsigned long -#else -#define UNSIGNED64 unsigned long long -#endif - -#ifndef U32 -#define U32 unsigned int -#endif - -/* exponent stored + 1024, hidden bit to left of decimal point */ -#define bias 1023 -#define bitstoright 52 -#define m1mask 0xf -#define hidden_bit 0x100000 -#ifdef LITTLE_ENDIAN_IEEE_DOUBLE -struct dblflt { - unsigned int m4: 16; - unsigned int m3: 16; - unsigned int m2: 16; - unsigned int m1: 4; - unsigned int e: 11; - unsigned int s: 1; -}; -#else -/* Big Endian IEEE Double Floats */ -struct dblflt { - unsigned int s: 1; - unsigned int e: 11; - unsigned int m1: 4; - unsigned int m2: 16; - unsigned int m3: 16; - unsigned int m4: 16; -}; -#endif -#define float_radix 2.147483648e9 - - -typedef UNSIGNED64 Bigit; -#define BIGSIZE 24 -#define MIN_E -1074 -#define MAX_FIVE 325 -#define B_P1 ((Bigit)1 << 52) - -typedef struct { - int l; - Bigit d[BIGSIZE]; -} Bignum; - -static Bignum R, S, MP, MM, five[MAX_FIVE]; -static Bignum S2, S3, S4, S5, S6, S7, S8, S9; -static int ruf, k, s_n, use_mp, qr_shift, sl, slr; - -static void mul10(Bignum *x); -static void big_short_mul(Bignum *x, Bigit y, Bignum *z); -/* -static void print_big(Bignum *x); -*/ -static int estimate(int n); -static void one_shift_left(int y, Bignum *z); -static void short_shift_left(Bigit x, int y, Bignum *z); -static void big_shift_left(Bignum *x, int y, Bignum *z); -static int big_comp(Bignum *x, Bignum *y); -static int sub_big(Bignum *x, Bignum *y, Bignum *z); -static void add_big(Bignum *x, Bignum *y, Bignum *z); -static int add_cmp(void); -static int qr(void); - -/*static int _PyFloat_Digits(char *buf, double v, int *signum);*/ -/*static void _PyFloat_DigitsInit(void);*/ - -#define ADD(x, y, z, k) {\ - Bigit x_add, z_add;\ - x_add = (x);\ - if ((k))\ - z_add = x_add + (y) + 1, (k) = (z_add <= x_add);\ - else\ - z_add = x_add + (y), (k) = (z_add < x_add);\ - (z) = z_add;\ -} - -#define SUB(x, y, z, b) {\ - Bigit x_sub, y_sub;\ - x_sub = (x); y_sub = (y);\ - if ((b))\ - (z) = x_sub - y_sub - 1, b = (y_sub >= x_sub);\ - else\ - (z) = x_sub - y_sub, b = (y_sub > x_sub);\ -} - -#define MUL(x, y, z, k) {\ - Bigit x_mul, low, high;\ - x_mul = (x);\ - low = (x_mul & 0xffffffff) * (y) + (k);\ - high = (x_mul >> 32) * (y) + (low >> 32);\ - (k) = high >> 32;\ - (z) = (low & 0xffffffff) | (high << 32);\ -} - -#define SLL(x, y, z, k) {\ - Bigit x_sll = (x);\ - (z) = (x_sll << (y)) | (k);\ - (k) = x_sll >> (64 - (y));\ -} - -static void -mul10(Bignum *x) -{ - int i, l; - Bigit *p, k; - - l = x->l; - for (i = l, p = &x->d[0], k = 0; i >= 0; i--) - MUL(*p, 10, *p++, k); - if (k != 0) - *p = k, x->l = l+1; -} - -static void -big_short_mul(Bignum *x, Bigit y, Bignum *z) -{ - int i, xl, zl; - Bigit *xp, *zp, k; - U32 high, low; - - xl = x->l; - xp = &x->d[0]; - zl = xl; - zp = &z->d[0]; - high = y >> 32; - low = y & 0xffffffff; - for (i = xl, k = 0; i >= 0; i--, xp++, zp++) { - Bigit xlow, xhigh, z0, t, c, z1; - xlow = *xp & 0xffffffff; - xhigh = *xp >> 32; - z0 = (xlow * low) + k; /* Cout is (z0 < k) */ - t = xhigh * low; - z1 = (xlow * high) + t; - c = (z1 < t); - t = z0 >> 32; - z1 += t; - c += (z1 < t); - *zp = (z1 << 32) | (z0 & 0xffffffff); - k = (xhigh * high) + (c << 32) + (z1 >> 32) + (z0 < k); - } - if (k != 0) - *zp = k, zl++; - z->l = zl; -} - -/* -static void -print_big(Bignum *x) -{ - int i; - Bigit *p; - - printf("#x"); - i = x->l; - p = &x->d[i]; - for (p = &x->d[i]; i >= 0; i--) { - Bigit b = *p--; - printf("%08x%08x", (int)(b >> 32), (int)(b & 0xffffffff)); - } -} -*/ - -static int -estimate(int n) -{ - if (n < 0) - return (int)(n*0.3010299956639812); - else - return 1+(int)(n*0.3010299956639811); -} - -static void -one_shift_left(int y, Bignum *z) -{ - int n, m, i; - Bigit *zp; - - n = y / 64; - m = y % 64; - zp = &z->d[0]; - for (i = n; i > 0; i--) *zp++ = 0; - *zp = (Bigit)1 << m; - z->l = n; -} - -static void -short_shift_left(Bigit x, int y, Bignum *z) -{ - int n, m, i, zl; - Bigit *zp; - - n = y / 64; - m = y % 64; - zl = n; - zp = &(z->d[0]); - for (i = n; i > 0; i--) *zp++ = 0; - if (m == 0) - *zp = x; - else { - Bigit high = x >> (64 - m); - *zp = x << m; - if (high != 0) - *++zp = high, zl++; - } - z->l = zl; -} - -static void -big_shift_left(Bignum *x, int y, Bignum *z) -{ - int n, m, i, xl, zl; - Bigit *xp, *zp, k; - - n = y / 64; - m = y % 64; - xl = x->l; - xp = &(x->d[0]); - zl = xl + n; - zp = &(z->d[0]); - for (i = n; i > 0; i--) *zp++ = 0; - if (m == 0) - for (i = xl; i >= 0; i--) *zp++ = *xp++; - else { - for (i = xl, k = 0; i >= 0; i--) - SLL(*xp++, m, *zp++, k); - if (k != 0) - *zp = k, zl++; - } - z->l = zl; -} - - -static int -big_comp(Bignum *x, Bignum *y) -{ - int i, xl, yl; - Bigit *xp, *yp; - - xl = x->l; - yl = y->l; - if (xl > yl) return 1; - if (xl < yl) return -1; - xp = &x->d[xl]; - yp = &y->d[xl]; - for (i = xl; i >= 0; i--, xp--, yp--) { - Bigit a = *xp; - Bigit b = *yp; - - if (a > b) return 1; - else if (a < b) return -1; - } - return 0; -} - -static int -sub_big(Bignum *x, Bignum *y, Bignum *z) -{ - int xl, yl, zl, b, i; - Bigit *xp, *yp, *zp; - - xl = x->l; - yl = y->l; - if (yl > xl) return 1; - xp = &x->d[0]; - yp = &y->d[0]; - zp = &z->d[0]; - - for (i = yl, b = 0; i >= 0; i--) - SUB(*xp++, *yp++, *zp++, b); - for (i = xl-yl; b && i > 0; i--) { - Bigit x_sub; - x_sub = *xp++; - *zp++ = x_sub - 1; - b = (x_sub == 0); - } - for (; i > 0; i--) *zp++ = *xp++; - if (b) return 1; - zl = xl; - while (*--zp == 0) zl--; - z->l = zl; - return 0; -} - -static void -add_big(Bignum *x, Bignum *y, Bignum *z) -{ - int xl, yl, k, i; - Bigit *xp, *yp, *zp; - - xl = x->l; - yl = y->l; - if (yl > xl) { - int tl; - Bignum *tn; - tl = xl; xl = yl; yl = tl; - tn = x; x = y; y = tn; - } - - xp = &x->d[0]; - yp = &y->d[0]; - zp = &z->d[0]; - - for (i = yl, k = 0; i >= 0; i--) - ADD(*xp++, *yp++, *zp++, k); - for (i = xl-yl; k && i > 0; i--) { - Bigit z_add; - z_add = *xp++ + 1; - k = (z_add == 0); - *zp++ = z_add; - } - for (; i > 0; i--) *zp++ = *xp++; - if (k) - *zp = 1, z->l = xl+1; - else - z->l = xl; -} - -static int -add_cmp() -{ - int rl, ml, sl, suml; - static Bignum sum; - - rl = R.l; - ml = (use_mp ? MP.l : MM.l); - sl = S.l; - - suml = rl >= ml ? rl : ml; - if ((sl > suml+1) || ((sl == suml+1) && (S.d[sl] > 1))) return -1; - if (sl < suml) return 1; - - add_big(&R, (use_mp ? &MP : &MM), &sum); - return big_comp(&sum, &S); -} - -static int -qr() -{ - if (big_comp(&R, &S5) < 0) - if (big_comp(&R, &S2) < 0) - if (big_comp(&R, &S) < 0) - return 0; - else { - sub_big(&R, &S, &R); - return 1; - } - else if (big_comp(&R, &S3) < 0) { - sub_big(&R, &S2, &R); - return 2; - } - else if (big_comp(&R, &S4) < 0) { - sub_big(&R, &S3, &R); - return 3; - } - else { - sub_big(&R, &S4, &R); - return 4; - } - else if (big_comp(&R, &S7) < 0) - if (big_comp(&R, &S6) < 0) { - sub_big(&R, &S5, &R); - return 5; - } - else { - sub_big(&R, &S6, &R); - return 6; - } - else if (big_comp(&R, &S9) < 0) - if (big_comp(&R, &S8) < 0) { - sub_big(&R, &S7, &R); - return 7; - } - else { - sub_big(&R, &S8, &R); - return 8; - } - else { - sub_big(&R, &S9, &R); - return 9; - } -} - -#define OUTDIG(d) { *buf++ = (d) + '0'; *buf = 0; return k; } - -int -_PyFloat_Digits(char *buf, double v, int *signum) -{ - struct dblflt *x; - int sign, e, f_n, m_n, i, d, tc1, tc2; - Bigit f; - - /* decompose float into sign, mantissa & exponent */ - x = (struct dblflt *)&v; - sign = x->s; - e = x->e; - f = (Bigit)(x->m1 << 16 | x->m2) << 32 | (U32)(x->m3 << 16 | x->m4); - if (e != 0) { - e = e - bias - bitstoright; - f |= (Bigit)hidden_bit << 32; - } - else if (f != 0) - /* denormalized */ - e = 1 - bias - bitstoright; - - *signum = sign; - if (f == 0) { - *buf++ = '0'; - *buf = 0; - return 0; - } - - ruf = !(f & 1); /* ruf = (even? f) */ - - /* Compute the scaling factor estimate, k */ - if (e > MIN_E) - k = estimate(e+52); - else { - int n; - Bigit y; - - for (n = e+52, y = (Bigit)1 << 52; f < y; n--) y >>= 1; - k = estimate(n); - } - - if (e >= 0) - if (f != B_P1) - use_mp = 0, f_n = e+1, s_n = 1, m_n = e; - else - use_mp = 1, f_n = e+2, s_n = 2, m_n = e; - else - if ((e == MIN_E) || (f != B_P1)) - use_mp = 0, f_n = 1, s_n = 1-e, m_n = 0; - else - use_mp = 1, f_n = 2, s_n = 2-e, m_n = 0; - - /* Scale it! */ - if (k == 0) { - short_shift_left(f, f_n, &R); - one_shift_left(s_n, &S); - one_shift_left(m_n, &MM); - if (use_mp) one_shift_left(m_n+1, &MP); - qr_shift = 1; - } - else if (k > 0) { - s_n += k; - if (m_n >= s_n) - f_n -= s_n, m_n -= s_n, s_n = 0; - else - f_n -= m_n, s_n -= m_n, m_n = 0; - short_shift_left(f, f_n, &R); - big_shift_left(&five[k-1], s_n, &S); - one_shift_left(m_n, &MM); - if (use_mp) one_shift_left(m_n+1, &MP); - qr_shift = 0; - } - else { - Bignum *power = &five[-k-1]; - - s_n += k; - big_short_mul(power, f, &S); - big_shift_left(&S, f_n, &R); - one_shift_left(s_n, &S); - big_shift_left(power, m_n, &MM); - if (use_mp) big_shift_left(power, m_n+1, &MP); - qr_shift = 1; - } - - /* fixup */ - if (add_cmp() <= -ruf) { - k--; - mul10(&R); - mul10(&MM); - if (use_mp) mul10(&MP); - } - - /* - printf("\nk = %d\n", k); - printf("R = "); print_big(&R); - printf("\nS = "); print_big(&S); - printf("\nM- = "); print_big(&MM); - if (use_mp) printf("\nM+ = "), print_big(&MP); - putchar('\n'); - fflush(0); - */ - - if (qr_shift) { - sl = s_n / 64; - slr = s_n % 64; - } - else { - big_shift_left(&S, 1, &S2); - add_big(&S2, &S, &S3); - big_shift_left(&S2, 1, &S4); - add_big(&S4, &S, &S5); - add_big(&S4, &S2, &S6); - add_big(&S4, &S3, &S7); - big_shift_left(&S4, 1, &S8); - add_big(&S8, &S, &S9); - } - -again: - if (qr_shift) { /* Take advantage of the fact that S = (ash 1 s_n) */ - if (R.l < sl) - d = 0; - else if (R.l == sl) { - Bigit *p; - - p = &R.d[sl]; - d = *p >> slr; - *p &= ((Bigit)1 << slr) - 1; - for (i = sl; (i > 0) && (*p == 0); i--) p--; - R.l = i; - } - else { - Bigit *p; - - p = &R.d[sl+1]; - d = *p << (64 - slr) | *(p-1) >> slr; - p--; - *p &= ((Bigit)1 << slr) - 1; - for (i = sl; (i > 0) && (*p == 0); i--) p--; - R.l = i; - } - } - else /* We need to do quotient-remainder */ - d = qr(); - - tc1 = big_comp(&R, &MM) < ruf; - tc2 = add_cmp() > -ruf; - if (!tc1) - if (!tc2) { - mul10(&R); - mul10(&MM); - if (use_mp) mul10(&MP); - *buf++ = d + '0'; - goto again; - } - else - OUTDIG(d+1) - else - if (!tc2) - OUTDIG(d) - else { - big_shift_left(&R, 1, &MM); - if (big_comp(&MM, &S) == -1) - OUTDIG(d) - else - OUTDIG(d+1) - } -} - -void -_PyFloat_DigitsInit() -{ - int n, i, l; - Bignum *b; - Bigit *xp, *zp, k; - - five[0].l = l = 0; - five[0].d[0] = 5; - for (n = MAX_FIVE-1, b = &five[0]; n > 0; n--) { - xp = &b->d[0]; - b++; - zp = &b->d[0]; - for (i = l, k = 0; i >= 0; i--) - MUL(*xp++, 5, *zp++, k); - if (k != 0) - *zp = k, l++; - b->l = l; - } - - /* - for (n = 1, b = &five[0]; n <= MAX_FIVE; n++) { - big_shift_left(b++, n, &R); - print_big(&R); - putchar('\n'); - } - fflush(0); - */ -} Index: Objects/intobject.c =================================================================== --- Objects/intobject.c (Revision 61828) +++ Objects/intobject.c (Arbeitskopie) @@ -1143,9 +1143,17 @@ return NULL; } +static PyObject * +int_is_finite(PyObject *v) +{ + Py_RETURN_TRUE; +} + static PyMethodDef int_methods[] = { {"conjugate", (PyCFunction)int_int, METH_NOARGS, "Returns self, the complex conjugate of any int."}, + {"is_finite", (PyCFunction)int_is_finite, METH_NOARGS, + "Returns always True."}, {"__trunc__", (PyCFunction)int_int, METH_NOARGS, "Truncating an Integral returns itself."}, {"__getnewargs__", (PyCFunction)int_getnewargs, METH_NOARGS}, Index: Objects/longobject.c =================================================================== --- Objects/longobject.c (Revision 61828) +++ Objects/longobject.c (Arbeitskopie) @@ -3428,9 +3428,17 @@ return NULL; } +static PyObject * +long_is_finite(PyObject *v) +{ + Py_RETURN_TRUE; +} + static PyMethodDef long_methods[] = { {"conjugate", (PyCFunction)long_long, METH_NOARGS, "Returns self, the complex conjugate of any long."}, + {"is_finite", (PyCFunction)long_is_finite, METH_NOARGS, + "Returns always True."}, {"__trunc__", (PyCFunction)long_long, METH_NOARGS, "Truncating an Integral returns itself."}, {"__getnewargs__", (PyCFunction)long_getnewargs, METH_NOARGS}, Index: Objects/floatobject.c =================================================================== --- Objects/floatobject.c (Revision 61828) +++ Objects/floatobject.c (Arbeitskopie) @@ -12,10 +12,12 @@ #include "formatter_string.h" -#if !defined(__STDC__) -extern double fmod(double, double); -extern double pow(double, double); +/* global NAN and INF objects */ +#ifdef Py_NAN +PyObject *PyFloat_NAN = NULL; #endif +PyObject *PyFloat_PINF = NULL; +PyObject *PyFloat_NINF = NULL; #ifdef _OSF_SOURCE /* OSF1 5.1 doesn't make this available with XOPEN_SOURCE_EXTENDED defined */ @@ -247,11 +249,11 @@ p++; } if (PyOS_strnicmp(p, "inf", 4) == 0) { - return PyFloat_FromDouble(sign * Py_HUGE_VAL); + Py_RETURN_INF(sign); } #ifdef Py_NAN if(PyOS_strnicmp(p, "nan", 4) == 0) { - return PyFloat_FromDouble(Py_NAN); + Py_RETURN_NAN; } #endif PyOS_snprintf(buffer, sizeof(buffer), @@ -401,110 +403,6 @@ format_float(buf, 100, v, precision); } -#ifdef Py_BROKEN_REPR -/* The following function is based on Tcl_PrintDouble, - * from tclUtil.c. - */ - -#define is_infinite(d) ( (d) > DBL_MAX || (d) < -DBL_MAX ) -#define is_nan(d) ((d) != (d)) - -static void -format_double_repr(char *dst, double value) -{ - char *p, c; - int exp; - int signum; - char buffer[30]; - - /* - * Handle NaN. - */ - - if (is_nan(value)) { - strcpy(dst, "nan"); - return; - } - - /* - * Handle infinities. - */ - - if (is_infinite(value)) { - if (value < 0) { - strcpy(dst, "-inf"); - } else { - strcpy(dst, "inf"); - } - return; - } - - /* - * Ordinary (normal and denormal) values. - */ - - exp = _PyFloat_Digits(buffer, value, &signum)+1; - if (signum) { - *dst++ = '-'; - } - p = buffer; - if (exp < -3 || exp > 17) { - /* - * E format for numbers < 1e-3 or >= 1e17. - */ - - *dst++ = *p++; - c = *p; - if (c != '\0') { - *dst++ = '.'; - while (c != '\0') { - *dst++ = c; - c = *++p; - } - } - sprintf(dst, "e%+d", exp-1); - } else { - /* - * F format for others. - */ - - if (exp <= 0) { - *dst++ = '0'; - } - c = *p; - while (exp-- > 0) { - if (c != '\0') { - *dst++ = c; - c = *++p; - } else { - *dst++ = '0'; - } - } - *dst++ = '.'; - if (c == '\0') { - *dst++ = '0'; - } else { - while (++exp < 0) { - *dst++ = '0'; - } - while (c != '\0') { - *dst++ = c; - c = *++p; - } - } - *dst++ = '\0'; - } -} - -static void -format_float_repr(char *buf, PyFloatObject *v) -{ - assert(PyFloat_Check(v)); - format_double_repr(buf, PyFloat_AS_DOUBLE(v)); -} - -#endif /* Py_BROKEN_REPR */ - /* Macro and helper that convert PyObject obj to a C double and store the value in dbl; this replaces the functionality of the coercion slot function. If conversion to double raises an exception, obj is @@ -589,13 +487,8 @@ static PyObject * float_repr(PyFloatObject *v) { -#ifdef Py_BROKEN_REPR - char buf[30]; - format_float_repr(buf, v); -#else char buf[100]; format_float(buf, sizeof(buf), v, PREC_REPR); -#endif return PyString_FromString(buf); } @@ -885,10 +778,13 @@ double a,b; CONVERT_TO_DOUBLE(v, a); CONVERT_TO_DOUBLE(w, b); +#ifdef Py_NAN if (b == 0.0) { - PyErr_SetString(PyExc_ZeroDivisionError, "float division"); + PyErr_SetString(PyExc_ZeroDivisionError, + "float division"); return NULL; } +#endif PyFPE_START_PROTECT("divide", return 0) a = a / b; PyFPE_END_PROTECT(a) @@ -904,10 +800,13 @@ if (Py_DivisionWarningFlag >= 2 && PyErr_Warn(PyExc_DeprecationWarning, "classic float division") < 0) return NULL; +#ifdef Py_NAN if (b == 0.0) { - PyErr_SetString(PyExc_ZeroDivisionError, "float division"); + PyErr_SetString(PyExc_ZeroDivisionError, + "float division"); return NULL; } +#endif PyFPE_START_PROTECT("divide", return 0) a = a / b; PyFPE_END_PROTECT(a) @@ -919,12 +818,15 @@ { double vx, wx; double mod; - CONVERT_TO_DOUBLE(v, vx); - CONVERT_TO_DOUBLE(w, wx); + CONVERT_TO_DOUBLE(v, vx); + CONVERT_TO_DOUBLE(w, wx); +#ifdef Py_NAN if (wx == 0.0) { - PyErr_SetString(PyExc_ZeroDivisionError, "float modulo"); + PyErr_SetString(PyExc_ZeroDivisionError, + "float modulo"); return NULL; } +#endif PyFPE_START_PROTECT("modulo", return 0) mod = fmod(vx, wx); /* note: checking mod*wx < 0 is incorrect -- underflows to @@ -1028,6 +930,9 @@ } return PyFloat_FromDouble(0.0); } + if (iv == 1.0) { /* 1**w is 1, even 1**inf and 1**nan */ + return PyFloat_FromDouble(1.0); + } if (iv < 0.0) { /* Whether this is an error is a mess, and bumps into libm * bugs so we have to figure it out ourselves. @@ -1119,6 +1024,55 @@ } static PyObject * +float_is_integer(PyObject *v) +{ + double x = PyFloat_AsDouble(v); + PyObject *o; + + if (x == -1.0 && PyErr_Occurred()) + return NULL; + if (!Py_IS_FINITE(x)) + Py_RETURN_FALSE; + PyFPE_START_PROTECT("is_integer", return NULL) + o = (floor(x) == x) ? Py_True : Py_False; + PyFPE_END_PROTECT(x) + if (errno != 0) { + PyErr_SetFromErrno(errno == ERANGE ? PyExc_OverflowError : + PyExc_ValueError); + return NULL; + } + Py_INCREF(o); + return o; +} + +static PyObject * +float_is_inf(PyObject *v) +{ + double x = PyFloat_AsDouble(v); + if (x == -1.0 && PyErr_Occurred()) + return NULL; + return PyBool_FromLong((long)Py_IS_INFINITY(x)); +} + +static PyObject * +float_is_nan(PyObject *v) +{ + double x = PyFloat_AsDouble(v); + if (x == -1.0 && PyErr_Occurred()) + return NULL; + return PyBool_FromLong((long)Py_IS_NAN(x)); +} + +static PyObject * +float_is_finite(PyObject *v) +{ + double x = PyFloat_AsDouble(v); + if (x == -1.0 && PyErr_Occurred()) + return NULL; + return PyBool_FromLong((long)Py_IS_FINITE(x)); +} + +static PyObject * float_trunc(PyObject *v) { double x = PyFloat_AsDouble(v); @@ -1476,12 +1430,20 @@ static PyMethodDef float_methods[] = { - {"conjugate", (PyCFunction)float_float, METH_NOARGS, + {"conjugate", (PyCFunction)float_float, METH_NOARGS, "Returns self, the complex conjugate of any float."}, {"__trunc__", (PyCFunction)float_trunc, METH_NOARGS, "Returns the Integral closest to x between 0 and x."}, {"as_integer_ratio", (PyCFunction)float_as_integer_ratio, METH_NOARGS, float_as_integer_ratio_doc}, + {"is_integer", (PyCFunction)float_is_integer, METH_NOARGS, + "Returns True if the float is an integer."}, + {"is_inf", (PyCFunction)float_is_inf, METH_NOARGS, + "Returns True if the float is positive or negative infinite."}, + {"is_finite", (PyCFunction)float_is_finite, METH_NOARGS, + "Returns True if the float is finite, neither infinite nor NaN."}, + {"is_nan", (PyCFunction)float_is_nan, METH_NOARGS, + "Returns True if the float is not a number (NaN)."}, {"__getnewargs__", (PyCFunction)float_getnewargs, METH_NOARGS}, {"__getformat__", (PyCFunction)float_getformat, METH_O|METH_CLASS, float_getformat_doc}, @@ -1642,13 +1604,29 @@ double_format = detected_double_format; float_format = detected_float_format; -#ifdef Py_BROKEN_REPR - /* Initialize floating point repr */ - _PyFloat_DigitsInit(); + /* Init float info */ + if (FloatInfoType.tp_name == 0) + PyStructSequence_InitType(&FloatInfoType, &floatinfo_desc); + + /* static floats */ +#define static_float(var, value) \ + if (var == NULL) \ + var = PyFloat_FromDouble(value); \ + else { \ + assert(PyFloat_CheckExact(var)); \ + _Py_NewReference(var); \ + } + +#ifdef Py_NAN + static_float(PyFloat_NAN, Py_NAN); #endif /* Init float info */ if (FloatInfoType.tp_name == 0) PyStructSequence_InitType(&FloatInfoType, &floatinfo_desc); + static_float(PyFloat_PINF, Py_HUGE_VAL); + static_float(PyFloat_NINF, -Py_HUGE_VAL); + +#undef static_float } void @@ -1709,6 +1687,12 @@ size_t bc, bf; /* block count, number of freed blocks */ size_t fsum; /* total unfreed floats per block */ +#ifdef Py_NAN + Py_XDECREF(PyFloat_NAN); +#endif + Py_XDECREF(PyFloat_PINF); + Py_XDECREF(PyFloat_NINF); + PyFloat_CompactFreeList(&bc, &bf, &fsum); if (!Py_VerboseFlag) Index: Misc/NEWS =================================================================== --- Misc/NEWS (Revision 61828) +++ Misc/NEWS (Arbeitskopie) @@ -11,7 +11,17 @@ Core and builtins ----------------- - + +- Added phase(z) -> phi, polar(z) -> r, phi and rect(r, phi) -> z to the cmath + module. + +- Issue #XXXX: Four new methods were added to the math and cmath modules: + acosh, asinh, atanh and log1p. Replacement implementations for platforms + without the four functions and copysign in libm were added to a new file + Python/pymath.c. + +- Added methods is_inf, is_nan and is_integer to float type. + - Issue #1477: With narrow Unicode builds, the unicode escape sequence \Uxxxxxxxx did not accept values outside the Basic Multilingual Plane. This affected raw unicode literals and the 'raw-unicode-escape' codec. Now Index: PC/pyconfig.h =================================================================== --- PC/pyconfig.h (Revision 61828) +++ PC/pyconfig.h (Arbeitskopie) @@ -207,12 +207,13 @@ #endif /* MS_WIN32 && !MS_WIN64 */ typedef int pid_t; -#define hypot _hypot #include #define Py_IS_NAN _isnan #define Py_IS_INFINITY(X) (!_finite(X) && !_isnan(X)) #define Py_IS_FINITE(X) _finite(X) +#define copysign _copysign +#define hypot _hypot #endif /* _MSC_VER */ @@ -392,7 +393,7 @@ /* Fairly standard from here! */ /* Define to 1 if you have the `copysign' function. */ -/* #define HAVE_COPYSIGN 1*/ +#define HAVE_COPYSIGN 1 /* Define to 1 if you have the `isinf' function. */ #define HAVE_ISINF 1 Index: Doc/library/cmath.rst =================================================================== --- Doc/library/cmath.rst (Revision 61828) +++ Doc/library/cmath.rst (Arbeitskopie) @@ -14,9 +14,82 @@ floating-point number, respectively, and the function is then applied to the result of the conversion. -The functions are: +.. 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. + +Complex coordinates +------------------- + +Complex numbers can be expressed by two important coordinate systems. +Python's :class:`complex` type uses rectangular coordinates where a number +on the complex plain is defined by two floats, the real part and the imaginary +part. + +Definition:: + + z = x + 1j * y + + x := real(z) + y := imag(z) + +In engineering the polar coordinate system is popular for complex numbers. In +polar coordinates a complex number is defined by the radius *r* and the phase +angle *φ*. The radius *r* is the absolute value of the complex, which can be +viewed as distance from (0, 0). The radius *r* is always 0 or a positive float. +The phase angle *φ* is the counter clockwise angle from the positive x axis, +e.g. *1* has the angle *0*, *1j* has the angle *π/2* and *-1* the angle *-π*. + +.. note:: + While :func:`phase` and func:`polar` return *+π* for a negative real they + may return *-π* for a complex with a very small negative imaginary + part, e.g. *-1-1E-300j*. + + +Definition:: + + z = r * exp(1j * φ) + z = r * cis(φ) + + r := abs(z) := sqrt(real(z)**2 + imag(z)**2) + phi := phase(z) := atan2(imag(z), real(z)) + cis(φ) := cos(φ) + 1j * sin(φ) + + +.. function:: phase(x) + + Return phase, also known as the argument, of a complex. + + .. versionadded:: 2.6 + + +.. function:: polar(x) + + Convert a :class:`complex` from rectangular coordinates to polar + coordinates. The function returns a tuple with the two elements + *r* and *phi*. *r* is the distance from 0 and *phi* the phase + angle. + + .. versionadded:: 2.6 + + +.. function:: rect(r, phi) + + Convert from polar coordinates to rectangular coordinates and return + a :class:`complex`. + + .. versionadded:: 2.6 + + + +cmath functions +--------------- + .. function:: acos(x) Return the arc cosine of *x*. There are two branch cuts: One extends right from @@ -37,32 +110,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*. @@ -78,6 +156,21 @@ Return the exponential value ``e**x``. +.. function:: isinf(x) + + Return *True* if the real or the imaginary part of x is positive + or negative infinity. + + .. versionadded:: 2.6 + + +.. function:: isnan(x) + + Return *True* if the real or imaginary part of x is not a number (NaN). + + .. versionadded:: 2.6 + + .. function:: log(x[, base]) Returns the logarithm of *x* to the given *base*. If the *base* is not @@ -154,3 +247,4 @@ nothing's sign bit. In Iserles, A., and Powell, M. (eds.), The state of the art in numerical analysis. Clarendon Press (1987) pp165-211. + Index: Doc/library/math.rst =================================================================== --- Doc/library/math.rst (Revision 61828) +++ Doc/library/math.rst (Arbeitskopie) @@ -139,6 +139,14 @@ *base* argument added. +.. function:: log1p(x) + + Return the natural logarithm of *1+x* (base *e*). The + result is calculated in a way which is accurate for *x* near zero. + + .. versionadded:: 2.6 + + .. function:: log10(x) Return the base-10 logarithm of *x*. @@ -146,9 +154,13 @@ .. function:: pow(x, y) - Return ``x**y``. + Return ``x**y``. ``1.0**y`` returns *1.0*, even for ``1.0**nan``. ``0**y`` + returns *0.* for all positive *y*, *0* and *NAN*. + .. versionchanged:: 2.6 + The outcome of ``1**nan`` and ``0**nan`` was undefined. + .. function:: sqrt(x) Return the square root of *x*. @@ -197,6 +209,13 @@ Return the sine of *x* radians. +.. function:: asinh(x) + + Return the inverse hyperbolic sine of *x*, in radians. + + .. versionadded:: 2.6 + + .. function:: tan(x) Return the tangent of *x* radians. @@ -221,6 +240,13 @@ Return the hyperbolic cosine of *x*. +.. function:: acosh(x) + + Return the inverse hyperbolic cosine of *x*, in radians. + + .. versionadded:: 2.6 + + .. function:: sinh(x) Return the hyperbolic sine of *x*. @@ -230,6 +256,14 @@ Return the hyperbolic tangent of *x*. + +.. function:: atanh(x) + + Return the inverse hyperbolic tangent of *x*, in radians. + + .. versionadded:: 2.6 + + The module also defines two mathematical constants: @@ -242,6 +276,7 @@ The mathematical constant *e*. + .. note:: The :mod:`math` module consists mostly of thin wrappers around the platform C @@ -255,9 +290,17 @@ :exc:`OverflowError` isn't defined, and in cases where ``math.log(0)`` raises :exc:`OverflowError`, ``math.log(0L)`` may raise :exc:`ValueError` instead. + All functions return a quite *NaN* if at least one of the args is *NaN*. + Signaling *NaN*s raise an exception. The exception type still depends on the + platform and libm implementation. It's usually :exc:`ValueError` for *EDOM* + and :exc:`OverflowError` for errno *ERANGE*. + ..versionchanged:: 2.6 + In earlier versions of Python the outcome of an operation with NaN as + input depended on platform and libm implementation. + + .. seealso:: Module :mod:`cmath` Complex number versions of many of these functions. - Index: Lib/test/test_math.py =================================================================== --- Lib/test/test_math.py (Revision 61828) +++ Lib/test/test_math.py (Arbeitskopie) @@ -4,10 +4,46 @@ from test.test_support import run_unittest, verbose import unittest import math +import os +import sys -seps='1e-05' -eps = eval(seps) +eps = 1E-05 +NAN = float('nan') +INF = float('inf') +NINF = float('-inf') +# locate file with test values +if __name__ == '__main__': + file = sys.argv[0] +else: + file = __file__ +test_dir = os.path.dirname(file) or os.curdir +test_file = os.path.join(test_dir, 'cmath_testcases.txt') + +def parse_testfile(fname): + """Parse a file with test values + + Empty lines or lines starting with -- are ignored + yields id, fn, arg_real, arg_imag, exp_real, exp_imag + """ + with open(fname) as fp: + for line in fp: + # skip comment lines and blank lines + if line.startswith('--') or not line.strip(): + continue + + lhs, rhs = line.split('->') + id, fn, arg_real, arg_imag = lhs.split() + rhs_pieces = rhs.split() + exp_real, exp_imag = rhs_pieces[0], rhs_pieces[1] + flags = rhs_pieces[2:] + + yield (id, fn, + float(arg_real), float(arg_imag), + float(exp_real), float(exp_imag), + flags + ) + class MathTests(unittest.TestCase): def ftest(self, name, value, expected): @@ -28,19 +64,58 @@ self.ftest('acos(-1)', math.acos(-1), math.pi) self.ftest('acos(0)', math.acos(0), math.pi/2) self.ftest('acos(1)', math.acos(1), 0) + self.assertRaises(ValueError, math.acos, INF) + self.assertRaises(ValueError, math.acos, NINF) + self.assert_(math.isnan(math.acos(NAN))) + def testAcosh(self): + self.assertRaises(TypeError, math.acosh) + self.ftest('acosh(1)', math.acosh(1), 0) + self.ftest('acosh(2)', math.acosh(2), 1.3169578969248168) + self.assertRaises(ValueError, math.acosh, 0) + self.assertRaises(ValueError, math.acosh, -1) + self.assertEquals(math.acosh(INF), INF) + self.assertRaises(ValueError, math.acosh, NINF) + self.assert_(math.isnan(math.acosh(NAN))) + def testAsin(self): self.assertRaises(TypeError, math.asin) self.ftest('asin(-1)', math.asin(-1), -math.pi/2) self.ftest('asin(0)', math.asin(0), 0) self.ftest('asin(1)', math.asin(1), math.pi/2) + self.assertRaises(ValueError, math.asin, INF) + self.assertRaises(ValueError, math.asin, NINF) + self.assert_(math.isnan(math.asin(NAN))) + def testAsinh(self): + self.assertRaises(TypeError, math.asinh) + self.ftest('asinh(0)', math.asinh(0), 0) + self.ftest('asinh(1)', math.asinh(1), 0.88137358701954305) + self.ftest('asinh(-1)', math.asinh(-1), -0.88137358701954305) + self.assertEquals(math.asinh(INF), INF) + self.assertEquals(math.asinh(NINF), NINF) + self.assert_(math.isnan(math.asinh(NAN))) + def testAtan(self): self.assertRaises(TypeError, math.atan) self.ftest('atan(-1)', math.atan(-1), -math.pi/4) self.ftest('atan(0)', math.atan(0), 0) self.ftest('atan(1)', math.atan(1), math.pi/4) + self.ftest('atan(inf)', math.atan(INF), math.pi/2) + self.ftest('atan(-inf)', math.atan(-INF), -math.pi/2) + self.assert_(math.isnan(math.atan(NAN))) + def testAtanh(self): + self.assertRaises(TypeError, math.atan) + self.ftest('atanh(0)', math.atanh(0), 0) + self.ftest('atanh(0.5)', math.atanh(0.5), 0.54930614433405489) + self.ftest('atanh(-0.5)', math.atanh(-0.5), -0.54930614433405489) + self.assertRaises(ValueError, math.atanh, 1) + self.assertRaises(ValueError, math.atanh, -1) + self.assertRaises(ValueError, math.atanh, INF) + self.assertRaises(ValueError, math.atanh, NINF) + self.assert_(math.isnan(math.atanh(NAN))) + def testAtan2(self): self.assertRaises(TypeError, math.atan2) self.ftest('atan2(-1, 0)', math.atan2(-1, 0), -math.pi/2) @@ -61,6 +136,9 @@ self.ftest('ceil(-0.5)', math.ceil(-0.5), 0) self.ftest('ceil(-1.0)', math.ceil(-1.0), -1) self.ftest('ceil(-1.5)', math.ceil(-1.5), -1) + self.assertEquals(math.ceil(INF), INF) + self.assertEquals(math.ceil(NINF), NINF) + self.assert_(math.isnan(math.ceil(NAN))) class TestCeil(object): def __float__(self): @@ -75,17 +153,55 @@ self.assertRaises(TypeError, math.ceil, t) self.assertRaises(TypeError, math.ceil, t, 0) + if float.__getformat__("double").startswith("IEEE"): + def testCopysign(self): + self.assertRaises(TypeError, math.copysign) + # copysign should let us distinguish signs of zeros + self.assertEquals(copysign(1., 0.), 1.) + self.assertEquals(copysign(1., -0.), -1.) + self.assertEquals(copysign(INF, 0.), INF) + self.assertEquals(copysign(INF, -0.), NINF) + self.assertEquals(copysign(NINF, 0.), INF) + self.assertEquals(copysign(NINF, -0.), NINF) + # and of infinities + self.assertEquals(copysign(1., INF), 1.) + self.assertEquals(copysign(1., NINF), -1.) + self.assertEquals(copysign(INF, INF), INF) + self.assertEquals(copysign(INF, NINF), NINF) + self.assertEquals(copysign(NINF, INF), INF) + self.assertEquals(copysign(NINF, NINF), NINF) + self.assert_(math.isnan(copysign(NAN, 1.))) + self.assert_(math.isnan(copysign(NAN, INF))) + self.assert_(math.isnan(copysign(NAN, NINF))) + self.assert_(math.isnan(copysign(NAN, NAN))) + # copysign(INF, NAN) may be INF or it may be NINF, since + # we don't know whether the sign bit of NAN is set on any + # given platform. + self.assert_(math.isinf(copysign(INF, NAN))) + # similarly, copysign(2., NAN) could be 2. or -2. + self.assertEquals(abs(copysign(2., NAN)), 2.) + def testCos(self): self.assertRaises(TypeError, math.cos) self.ftest('cos(-pi/2)', math.cos(-math.pi/2), 0) self.ftest('cos(0)', math.cos(0), 1) self.ftest('cos(pi/2)', math.cos(math.pi/2), 0) self.ftest('cos(pi)', math.cos(math.pi), -1) + try: + self.assert_(math.isnan(math.cos(INF))) + self.assert_(math.isnan(math.cos(NINF))) + except ValueError: + self.assertRaises(ValueError, math.cos, INF) + self.assertRaises(ValueError, math.cos, NINF) + self.assert_(math.isnan(math.cos(NAN))) def testCosh(self): self.assertRaises(TypeError, math.cosh) self.ftest('cosh(0)', math.cosh(0), 1) self.ftest('cosh(2)-2*cosh(1)**2', math.cosh(2)-2*math.cosh(1)**2, -1) # Thanks to Lambert + self.assertEquals(math.cosh(INF), INF) + self.assertEquals(math.cosh(NINF), INF) + self.assert_(math.isnan(math.cosh(NAN))) def testDegrees(self): self.assertRaises(TypeError, math.degrees) @@ -98,6 +214,9 @@ self.ftest('exp(-1)', math.exp(-1), 1/math.e) self.ftest('exp(0)', math.exp(0), 1) self.ftest('exp(1)', math.exp(1), math.e) + self.assertEquals(math.exp(INF), INF) + self.assertEquals(math.exp(NINF), 0.) + self.assert_(math.isnan(math.exp(NAN))) def testFabs(self): self.assertRaises(TypeError, math.fabs) @@ -121,6 +240,9 @@ # This fails on some platforms - so check it here self.ftest('floor(1.23e167)', math.floor(1.23e167), 1.23e167) self.ftest('floor(-1.23e167)', math.floor(-1.23e167), -1.23e167) + self.assertEquals(math.ceil(INF), INF) + self.assertEquals(math.ceil(NINF), NINF) + self.assert_(math.isnan(math.floor(NAN))) class TestFloor(object): def __float__(self): @@ -143,6 +265,19 @@ self.ftest('fmod(-10,1)', math.fmod(-10,1), 0) self.ftest('fmod(-10,0.5)', math.fmod(-10,0.5), 0) self.ftest('fmod(-10,1.5)', math.fmod(-10,1.5), -1) + self.assert_(math.isnan(math.fmod(NAN, 1.))) + self.assert_(math.isnan(math.fmod(1., NAN))) + self.assert_(math.isnan(math.fmod(NAN, NAN))) + self.assertRaises(ValueError, math.fmod, 1., 0.) + self.assertRaises(ValueError, math.fmod, INF, 1.) + self.assertRaises(ValueError, math.fmod, NINF, 1.) + self.assertRaises(ValueError, math.fmod, INF, 0.) + self.assertEquals(math.fmod(3.0, INF), 3.0) + self.assertEquals(math.fmod(-3.0, INF), -3.0) + self.assertEquals(math.fmod(3.0, NINF), 3.0) + self.assertEquals(math.fmod(-3.0, NINF), -3.0) + self.assertEquals(math.fmod(0.0, 3.0), 0.0) + self.assertEquals(math.fmod(0.0, NINF), 0.0) def testFrexp(self): self.assertRaises(TypeError, math.frexp) @@ -157,10 +292,20 @@ testfrexp('frexp(1)', math.frexp(1), (0.5, 1)) testfrexp('frexp(2)', math.frexp(2), (0.5, 2)) + self.assertEquals(math.frexp(INF)[0], INF) + self.assertEquals(math.frexp(NINF)[0], NINF) + self.assert_(math.isnan(math.frexp(NAN)[0])) + def testHypot(self): self.assertRaises(TypeError, math.hypot) self.ftest('hypot(0,0)', math.hypot(0,0), 0) self.ftest('hypot(3,4)', math.hypot(3,4), 5) + self.assertEqual(math.hypot(NAN, INF), INF) + self.assertEqual(math.hypot(INF, NAN), INF) + self.assertEqual(math.hypot(NAN, NINF), INF) + self.assertEqual(math.hypot(NINF, NAN), INF) + self.assert_(math.isnan(math.hypot(1.0, NAN))) + self.assert_(math.isnan(math.hypot(NAN, -2.0))) def testLdexp(self): self.assertRaises(TypeError, math.ldexp) @@ -168,6 +313,13 @@ self.ftest('ldexp(1,1)', math.ldexp(1,1), 2) self.ftest('ldexp(1,-1)', math.ldexp(1,-1), 0.5) self.ftest('ldexp(-1,1)', math.ldexp(-1,1), -2) + self.assertRaises(OverflowError, math.ldexp, 1., 1000000) + self.assertRaises(OverflowError, math.ldexp, -1., 1000000) + self.assertEquals(math.ldexp(1., -1000000), 0.) + self.assertEquals(math.ldexp(-1., -1000000), -0.) + self.assertEquals(math.ldexp(INF, 30), INF) + self.assertEquals(math.ldexp(NINF, -213), NINF) + self.assert_(math.isnan(math.ldexp(NAN, 0))) def testLog(self): self.assertRaises(TypeError, math.log) @@ -177,12 +329,31 @@ self.ftest('log(32,2)', math.log(32,2), 5) self.ftest('log(10**40, 10)', math.log(10**40, 10), 40) self.ftest('log(10**40, 10**20)', math.log(10**40, 10**20), 2) + self.assertEquals(math.log(INF), INF) + self.assertRaises(ValueError, math.log, NINF) + self.assert_(math.isnan(math.log(NAN))) + def testLog1p(self): + self.assertRaises(TypeError, math.log1p) + self.ftest('log1p(1/e -1)', math.log1p(1/math.e-1), -1) + self.ftest('log1p(0)', math.log1p(0), 0) + self.ftest('log1p(e-1)', math.log1p(math.e-1), 1) + self.ftest('log1p(1)', math.log1p(1), math.log(2)) + self.assertEquals(math.log1p(INF), INF) + self.assertRaises(ValueError, math.log1p, NINF) + self.assert_(math.isnan(math.log1p(NAN))) + n= 2**90 + self.assertAlmostEquals(math.log1p(n), 62.383246250395075) + self.assertAlmostEquals(math.log1p(n), math.log1p(float(n))) + def testLog10(self): self.assertRaises(TypeError, math.log10) self.ftest('log10(0.1)', math.log10(0.1), -1) self.ftest('log10(1)', math.log10(1), 0) self.ftest('log10(10)', math.log10(10), 1) + self.assertEquals(math.log(INF), INF) + self.assertRaises(ValueError, math.log10, NINF) + self.assert_(math.isnan(math.log10(NAN))) def testModf(self): self.assertRaises(TypeError, math.modf) @@ -195,12 +366,35 @@ testmodf('modf(1.5)', math.modf(1.5), (0.5, 1.0)) testmodf('modf(-1.5)', math.modf(-1.5), (-0.5, -1.0)) + self.assertEquals(math.modf(INF), (0.0, INF)) + self.assertEquals(math.modf(NINF), (-0.0, NINF)) + + modf_nan = math.modf(NAN) + self.assert_(math.isnan(modf_nan[0])) + self.assert_(math.isnan(modf_nan[1])) + def testPow(self): self.assertRaises(TypeError, math.pow) self.ftest('pow(0,1)', math.pow(0,1), 0) self.ftest('pow(1,0)', math.pow(1,0), 1) self.ftest('pow(2,1)', math.pow(2,1), 2) self.ftest('pow(2,-1)', math.pow(2,-1), 0.5) + self.assertEqual(math.pow(INF, 1), INF) + self.assertEqual(math.pow(NINF, 1), NINF) + self.assertEqual((math.pow(1, INF)), 1.) + self.assertEqual((math.pow(1, NINF)), 1.) + self.assert_(math.isnan(math.pow(NAN, 1))) + self.assert_(math.isnan(math.pow(2, NAN))) + self.assert_(math.isnan(math.pow(0, NAN))) + self.assertEqual(math.pow(1, NAN), 1) + self.assertEqual(1**NAN, 1) + self.assertEqual(1**INF, 1) + self.assertEqual(1**NINF, 1) + self.assertEqual(1**0, 1) + self.assertEqual(1.**NAN, 1) + self.assertEqual(1.**INF, 1) + self.assertEqual(1.**NINF, 1) + self.assertEqual(1.**0, 1) def testRadians(self): self.assertRaises(TypeError, math.radians) @@ -213,29 +407,52 @@ self.ftest('sin(0)', math.sin(0), 0) self.ftest('sin(pi/2)', math.sin(math.pi/2), 1) self.ftest('sin(-pi/2)', math.sin(-math.pi/2), -1) + try: + self.assert_(math.isnan(math.sin(INF))) + self.assert_(math.isnan(math.sin(NINF))) + except ValueError: + self.assertRaises(ValueError, math.sin, INF) + self.assertRaises(ValueError, math.sin, NINF) + self.assert_(math.isnan(math.sin(NAN))) def testSinh(self): self.assertRaises(TypeError, math.sinh) self.ftest('sinh(0)', math.sinh(0), 0) self.ftest('sinh(1)**2-cosh(1)**2', math.sinh(1)**2-math.cosh(1)**2, -1) self.ftest('sinh(1)+sinh(-1)', math.sinh(1)+math.sinh(-1), 0) + self.assertEquals(math.sinh(INF), INF) + self.assertEquals(math.sinh(-INF), -INF) + self.assert_(math.isnan(math.sinh(NAN))) def testSqrt(self): self.assertRaises(TypeError, math.sqrt) self.ftest('sqrt(0)', math.sqrt(0), 0) self.ftest('sqrt(1)', math.sqrt(1), 1) self.ftest('sqrt(4)', math.sqrt(4), 2) + self.assertEquals(math.sqrt(INF), INF) + self.assertRaises(ValueError, math.sqrt, NINF) + self.assert_(math.isnan(math.sqrt(NAN))) def testTan(self): self.assertRaises(TypeError, math.tan) self.ftest('tan(0)', math.tan(0), 0) self.ftest('tan(pi/4)', math.tan(math.pi/4), 1) self.ftest('tan(-pi/4)', math.tan(-math.pi/4), -1) + try: + self.assert_(math.isnan(math.tan(INF))) + self.assert_(math.isnan(math.tan(NINF))) + except: + self.assertRaises(ValueError, math.tan, INF) + self.assertRaises(ValueError, math.tan, NINF) + self.assert_(math.isnan(math.tan(NAN))) def testTanh(self): self.assertRaises(TypeError, math.tanh) self.ftest('tanh(0)', math.tanh(0), 0) self.ftest('tanh(1)+tanh(-1)', math.tanh(1)+math.tanh(-1), 0) + self.ftest('tanh(inf)', math.tanh(INF), 1) + self.ftest('tanh(-inf)', math.tanh(NINF), -1) + self.assert_(math.isnan(math.tanh(NAN))) def test_trunc(self): self.assertEqual(math.trunc(1), 1) @@ -329,9 +546,27 @@ else: self.fail("sqrt(-1) didn't raise ValueError") + def test_testfile(self): + if not float.__getformat__("double").startswith("IEEE"): + return + for id, fn, ar, ai, er, ei, flags in parse_testfile(test_file): + # Skip if either the input or result is complex, or if + # flags is nonempty + if ai != 0. or ei != 0. or flags: + continue + if fn in ['rect', 'polar']: + # no real versions of rect, polar + continue + func = getattr(math, fn) + result = func(ar) + self.ftest("%s:%s(%r)" % (id, fn, ar), result, er) def test_main(): - run_unittest(MathTests) + from doctest import DocFileSuite + suite = unittest.TestSuite() + suite.addTest(unittest.makeSuite(MathTests)) + suite.addTest(DocFileSuite("ieee754.txt")) + run_unittest(suite) if __name__ == '__main__': test_main() Index: Lib/test/test_float.py =================================================================== --- Lib/test/test_float.py (Revision 61828) +++ Lib/test/test_float.py (Arbeitskopie) @@ -2,13 +2,13 @@ import unittest, struct import os from test import test_support +import math +from math import isinf, isnan +import operator -def isinf(x): - return x * 0.5 == x +INF = float("inf") +NAN = float("nan") -def isnan(x): - return x != x - class FormatFunctionsTestCase(unittest.TestCase): def setUp(self): @@ -206,6 +206,17 @@ self.assertEqual(str(1e300 * 1e300 * 0), "nan") self.assertEqual(str(-1e300 * 1e300 * 0), "nan") + def test_float_nan(self): + self.assert_(NAN.is_nan()) + self.failIf(INF.is_nan()) + self.failIf((0.).is_nan()) + + def test_float_inf(self): + self.assert_(INF.is_inf()) + self.failIf(NAN.is_inf()) + self.failIf((0.).is_inf()) + + def test_main(): test_support.run_unittest( FormatFunctionsTestCase, Index: Lib/test/test_cmath.py =================================================================== --- Lib/test/test_cmath.py (Revision 61828) +++ Lib/test/test_cmath.py (Arbeitskopie) @@ -1,7 +1,82 @@ from test.test_support import run_unittest +from test.test_math import parse_testfile, test_file import unittest +import os, sys import cmath, math +from cmath import phase, polar, rect, pi +INF = float('inf') +NAN = float('nan') + +complex_zeros = [complex(x, y) for x in [0.0, -0.0] for y in [0.0, -0.0]] +complex_infinities = [complex(x, y) for x, y in [ + (INF, 0.0), # 1st quadrant + (INF, 2.3), + (INF, INF), + (2.3, INF), + (0.0, INF), + (-0.0, INF), # 2nd quadrant + (-2.3, INF), + (-INF, INF), + (-INF, 2.3), + (-INF, 0.0), + (-INF, -0.0), # 3rd quadrant + (-INF, -2.3), + (-INF, -INF), + (-2.3, -INF), + (-0.0, -INF), + (0.0, -INF), # 4th quadrant + (2.3, -INF), + (INF, -INF), + (INF, -2.3), + (INF, -0.0) + ]] +complex_nans = [complex(x, y) for x, y in [ + (NAN, -INF), + (NAN, -2.3), + (NAN, -0.0), + (NAN, 0.0), + (NAN, 2.3), + (NAN, INF), + (-INF, NAN), + (-2.3, NAN), + (-0.0, NAN), + (0.0, NAN), + (2.3, NAN), + (INF, NAN) + ]] + +def almostEqualF(a, b, rel_err=2e-15, abs_err = 5e-323): + """Determine whether floating-point values a and b are equal to within + a (small) rounding error. The default values for rel_err and + abs_err are chosen to be suitable for platforms where a float is + represented by an IEEE 754 double. They allow an error of between + 9 and 19 ulps.""" + + # special values testing + if math.isnan(a): + return math.isnan(b) + if math.isinf(a): + return a == b + + # if both a and b are zero, check whether they have the same sign + # (in theory there are examples where it would be legitimate for a + # and b to have opposite signs; in practice these hardly ever + # occur). + if not a and not b: + return math.copysign(1., a) == math.copysign(1., b) + + # if a-b overflows, or b is infinite, return False. Again, in + # theory there are examples where a is within a few ulps of the + # max representable float, and then b could legitimately be + # infinite. In practice these examples are rare. + try: + absolute_error = abs(b-a) + except OverflowError: + return False + else: + return absolute_error <= max(abs_err, rel_err * abs(a)) + class CMathTests(unittest.TestCase): # list of all functions in cmath test_functions = [getattr(cmath, fname) for fname in [ @@ -12,24 +87,52 @@ test_functions.append(lambda x : cmath.log(x, 1729. + 0j)) test_functions.append(lambda x : cmath.log(14.-27j, x)) - def cAssertAlmostEqual(self, a, b, rel_eps = 1e-10, abs_eps = 1e-100): - """Check that two complex numbers are almost equal.""" - # the two complex numbers are considered almost equal if - # either the relative error is <= rel_eps or the absolute error - # is tiny, <= abs_eps. - if a == b == 0: - return - absolute_error = abs(a-b) - relative_error = absolute_error/max(abs(a), abs(b)) - if relative_error > rel_eps and absolute_error > abs_eps: - self.fail("%s and %s are not almost equal" % (a, b)) + def setUp(self): + self.test_values = open(test_file) + def tearDown(self): + self.test_values.close() + + def rAssertAlmostEqual(self, a, b, rel_err = 2e-15, abs_err = 5e-323): + """Check that two floating-point numbers are almost equal.""" + + # special values testing + if math.isnan(a): + if math.isnan(b): + return + self.fail("%s should be nan" % repr(b)) + + if math.isinf(a): + if a == b: + return + self.fail("finite result where infinity excpected: " + "expected %s, got %s" % (repr(a), repr(b))) + + if not a and not b: + if math.atan2(a, -1.) != math.atan2(b, -1.): + self.fail("zero has wrong sign: expected %s, got %s" % + (repr(a), repr(b))) + + # test passes if either the absolute error or the relative + # error is sufficiently small. The defaults amount to an + # error of between 9 ulps and 19 ulps on an IEEE-754 compliant + # machine. + + try: + absolute_error = abs(b-a) + except OverflowError: + pass + else: + if absolute_error <= max(abs_err, rel_err * abs(a)): + return + self.fail("%s and %s are not sufficiently close" % (repr(a), repr(b))) + def test_constants(self): e_expected = 2.71828182845904523536 pi_expected = 3.14159265358979323846 - self.assertAlmostEqual(cmath.pi, pi_expected, 9, + self.rAssertAlmostEqual(cmath.pi, pi_expected, 9, "cmath.pi is %s; should be %s" % (cmath.pi, pi_expected)) - self.assertAlmostEqual(cmath.e, e_expected, 9, + self.rAssertAlmostEqual(cmath.e, e_expected, 9, "cmath.e is %s; should be %s" % (cmath.e, e_expected)) def test_user_object(self): @@ -109,13 +212,13 @@ for f in self.test_functions: # usual usage - self.cAssertAlmostEqual(f(MyComplex(cx_arg)), f(cx_arg)) - self.cAssertAlmostEqual(f(MyComplexOS(cx_arg)), f(cx_arg)) + self.assertEqual(f(MyComplex(cx_arg)), f(cx_arg)) + self.assertEqual(f(MyComplexOS(cx_arg)), f(cx_arg)) # other combinations of __float__ and __complex__ - self.cAssertAlmostEqual(f(FloatAndComplex()), f(cx_arg)) - self.cAssertAlmostEqual(f(FloatAndComplexOS()), f(cx_arg)) - self.cAssertAlmostEqual(f(JustFloat()), f(flt_arg)) - self.cAssertAlmostEqual(f(JustFloatOS()), f(flt_arg)) + self.assertEqual(f(FloatAndComplex()), f(cx_arg)) + self.assertEqual(f(FloatAndComplexOS()), f(cx_arg)) + self.assertEqual(f(JustFloat()), f(flt_arg)) + self.assertEqual(f(JustFloatOS()), f(flt_arg)) # TypeError should be raised for classes not providing # either __complex__ or __float__, even if they provide # __int__, __long__ or __index__. An old-style class @@ -138,7 +241,7 @@ # functions, by virtue of providing a __float__ method for f in self.test_functions: for arg in [2, 2L, 2.]: - self.cAssertAlmostEqual(f(arg), f(arg.__float__())) + self.assertEqual(f(arg), f(arg.__float__())) # but strings should give a TypeError for f in self.test_functions: @@ -182,13 +285,202 @@ float_fn = getattr(math, fn) complex_fn = getattr(cmath, fn) for v in values: - self.cAssertAlmostEqual(float_fn(v), complex_fn(v)) + z = complex_fn(v) + self.rAssertAlmostEqual(float_fn(v), z.real) + self.assertEqual(0., z.imag) # test two-argument version of log with various bases for base in [0.5, 2., 10.]: for v in positive: - self.cAssertAlmostEqual(cmath.log(v, base), math.log(v, base)) + z = cmath.log(v, base) + self.rAssertAlmostEqual(math.log(v, base), z.real) + self.assertEqual(0., z.imag) + def test_specific_values(self): + if not float.__getformat__("double").startswith("IEEE"): + return + + def rect_complex(z): + """Wrapped version of rect that accepts a complex number instead of + two float arguments.""" + return cmath.rect(z.real, z.imag) + + def polar_complex(z): + """Wrapped version of polar that returns a complex number instead of + two floats.""" + return complex(*polar(z)) + + for id, fn, ar, ai, er, ei, flags in parse_testfile(test_file): + arg = complex(ar, ai) + expected = complex(er, ei) + if fn == 'rect': + function = rect_complex + elif fn == 'polar': + function = polar_complex + else: + function = getattr(cmath, fn) + if 'divide-by-zero' in flags or 'invalid' in flags: + try: + actual = function(arg) + except ValueError: + continue + else: + test_str = "%s: %s(complex(%r, %r))" % (id, fn, ar, ai) + self.fail('ValueError not raised in test %s' % test_str) + + if 'overflow' in flags: + try: + actual = function(arg) + except OverflowError: + continue + else: + test_str = "%s: %s(complex(%r, %r))" % (id, fn, ar, ai) + self.fail('OverflowError not raised in test %s' % test_str) + + actual = function(arg) + + if 'ignore-real-sign' in flags: + actual = complex(abs(actual.real), actual.imag) + expected = complex(abs(expected.real), expected.imag) + if 'ignore-imag-sign' in flags: + actual = complex(actual.real, abs(actual.imag)) + expected = complex(expected.real, abs(expected.imag)) + + # for the real part of the log function, we allow an + # absolute error of up to 2e-15. + if fn in ('log', 'log10'): + real_abs_err = 2e-15 + else: + real_abs_err = 5e-323 + + if not (almostEqualF(expected.real, actual.real, + abs_err = real_abs_err) and + almostEqualF(expected.imag, actual.imag)): + error_message = ( + "%s: %s(complex(%r, %r))\n" % (id, fn, ar, ai) + + "Expected: complex(%r, %r)\n" % + (expected.real, expected.imag) + + "Received: complex(%r, %r)\n" % + (actual.real, actual.imag) + + "Received value insufficiently close to expected value.") + self.fail(error_message) + + def assertCISEqual(self, a, b): + eps = 1E-7 + if abs(a[0] - b[0]) > eps or abs(a[1] - b[1]) > eps: + self.fail((a ,b)) + + def test_polar(self): + self.assertCISEqual(polar(0), (0., 0.)) + self.assertCISEqual(polar(1.), (1., 0.)) + self.assertCISEqual(polar(-1.), (1., pi)) + self.assertCISEqual(polar(1j), (1., pi/2)) + self.assertCISEqual(polar(-1j), (1., -pi/2)) + + def test_phase(self): + self.assertAlmostEqual(phase(0), 0.) + self.assertAlmostEqual(phase(1.), 0.) + self.assertAlmostEqual(phase(-1.), pi) + self.assertAlmostEqual(phase(-1.+1E-300j), pi) + self.assertAlmostEqual(phase(-1.-1E-300j), -pi) + self.assertAlmostEqual(phase(1j), pi/2) + self.assertAlmostEqual(phase(-1j), -pi/2) + + # zeros + self.assertEqual(phase(complex(0.0, 0.0)), 0.0) + self.assertEqual(phase(complex(0.0, -0.0)), -0.0) + self.assertEqual(phase(complex(-0.0, 0.0)), pi) + self.assertEqual(phase(complex(-0.0, -0.0)), -pi) + + # infinities + self.assertAlmostEqual(phase(complex(-INF, -0.0)), -pi) + self.assertAlmostEqual(phase(complex(-INF, -2.3)), -pi) + self.assertAlmostEqual(phase(complex(-INF, -INF)), -0.75*pi) + self.assertAlmostEqual(phase(complex(-2.3, -INF)), -pi/2) + self.assertAlmostEqual(phase(complex(-0.0, -INF)), -pi/2) + self.assertAlmostEqual(phase(complex(0.0, -INF)), -pi/2) + self.assertAlmostEqual(phase(complex(2.3, -INF)), -pi/2) + self.assertAlmostEqual(phase(complex(INF, -INF)), -pi/4) + self.assertEqual(phase(complex(INF, -2.3)), -0.0) + self.assertEqual(phase(complex(INF, -0.0)), -0.0) + self.assertEqual(phase(complex(INF, 0.0)), 0.0) + self.assertEqual(phase(complex(INF, 2.3)), 0.0) + self.assertAlmostEqual(phase(complex(INF, INF)), pi/4) + self.assertAlmostEqual(phase(complex(2.3, INF)), pi/2) + self.assertAlmostEqual(phase(complex(0.0, INF)), pi/2) + self.assertAlmostEqual(phase(complex(-0.0, INF)), pi/2) + self.assertAlmostEqual(phase(complex(-2.3, INF)), pi/2) + self.assertAlmostEqual(phase(complex(-INF, INF)), 0.75*pi) + self.assertAlmostEqual(phase(complex(-INF, 2.3)), pi) + self.assertAlmostEqual(phase(complex(-INF, 0.0)), pi) + + # real or imaginary part NaN + for z in complex_nans: + self.assert_(math.isnan(phase(z))) + + def test_abs(self): + # zeros + for z in complex_zeros: + self.assertEqual(abs(z), 0.0) + + # infinities + for z in complex_infinities: + self.assertEqual(abs(z), INF) + + # real or imaginary part NaN + self.assertEqual(abs(complex(NAN, -INF)), INF) + self.assert_(math.isnan(abs(complex(NAN, -2.3)))) + self.assert_(math.isnan(abs(complex(NAN, -0.0)))) + self.assert_(math.isnan(abs(complex(NAN, 0.0)))) + self.assert_(math.isnan(abs(complex(NAN, 2.3)))) + self.assertEqual(abs(complex(NAN, INF)), INF) + self.assertEqual(abs(complex(-INF, NAN)), INF) + self.assert_(math.isnan(abs(complex(-2.3, NAN)))) + self.assert_(math.isnan(abs(complex(-0.0, NAN)))) + self.assert_(math.isnan(abs(complex(0.0, NAN)))) + self.assert_(math.isnan(abs(complex(2.3, NAN)))) + self.assertEqual(abs(complex(INF, NAN)), INF) + self.assert_(math.isnan(abs(complex(NAN, NAN)))) + + # result overflows + if float.__getformat__("double").startswith("IEEE"): + self.assertRaises(OverflowError, abs, complex(1.4e308, 1.4e308)) + + def assertCEqual(self, a, b): + eps = 1E-7 + if abs(a.real - b[0]) > eps or abs(a.imag - b[1]) > eps: + self.fail((a ,b)) + + def test_rect(self): + self.assertCEqual(rect(0, 0), (0, 0)) + self.assertCEqual(rect(1, 0), (1., 0)) + self.assertCEqual(rect(1, -pi), (-1., 0)) + self.assertCEqual(rect(1, pi/2), (0, 1.)) + self.assertCEqual(rect(1, -pi/2), (0, -1.)) + + def test_isnan(self): + self.failIf(cmath.isnan(1)) + self.failIf(cmath.isnan(1j)) + self.failIf(cmath.isnan(INF)) + self.assert_(cmath.isnan(NAN)) + self.assert_(cmath.isnan(complex(NAN, 0))) + self.assert_(cmath.isnan(complex(0, NAN))) + self.assert_(cmath.isnan(complex(NAN, NAN))) + self.assert_(cmath.isnan(complex(NAN, INF))) + self.assert_(cmath.isnan(complex(INF, NAN))) + + def test_isinf(self): + self.failIf(cmath.isinf(1)) + self.failIf(cmath.isinf(1j)) + self.failIf(cmath.isinf(NAN)) + self.assert_(cmath.isinf(INF)) + self.assert_(cmath.isinf(complex(INF, 0))) + self.assert_(cmath.isinf(complex(0, INF))) + self.assert_(cmath.isinf(complex(INF, INF))) + self.assert_(cmath.isinf(complex(NAN, INF))) + self.assert_(cmath.isinf(complex(INF, NAN))) + + def test_main(): run_unittest(CMathTests) Index: Makefile.pre.in =================================================================== --- Makefile.pre.in (Revision 61828) +++ Makefile.pre.in (Arbeitskopie) @@ -272,6 +272,7 @@ Python/peephole.o \ Python/pyarena.o \ Python/pyfpe.o \ + Python/pymath.o \ Python/pystate.o \ Python/pythonrun.o \ Python/structmember.o \ @@ -594,6 +595,7 @@ Include/pydebug.h \ Include/pyerrors.h \ Include/pyfpe.h \ + Include/pymath.h \ Include/pygetopt.h \ Include/pymem.h \ Include/pyport.h \ Index: Modules/cmathmodule.c =================================================================== --- Modules/cmathmodule.c (Revision 61828) +++ Modules/cmathmodule.c (Arbeitskopie) @@ -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 -#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 -/* 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_LN10 +#define M_LN10 (2.302585092994045684) /* natural log of 10 */ +#endif +/* + 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 @@ "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 @@ 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 @@ "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,21 +284,100 @@ 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, "atan(x)\n" "\n" "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 @@ 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 @@ "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 @@ "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 @@ "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,38 +638,82 @@ "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" "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 @@ "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 @@ 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 @@ "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 @@ "\n" "Return the hyperbolic tangent of x."); + static PyObject * cmath_log(PyObject *self, PyObject *args) { @@ -325,7 +923,6 @@ 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 @@ 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,7 +989,152 @@ 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" "functions for complex numbers."); @@ -401,8 +1149,13 @@ {"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 @@ 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)); } Index: Modules/mathmodule.c =================================================================== --- Modules/mathmodule.c (Revision 61828) +++ Modules/mathmodule.c (Arbeitskopie) @@ -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,28 +95,97 @@ 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(PyObject *arg, double (*func) (double)) +math_1(PyObject *arg, 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 PyFloat_FromDouble(x); + return PyFloat_FromDouble(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_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); @@ -81,19 +193,30 @@ 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); @@ -103,55 +226,49 @@ }\ 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(ceil, ceil, +FUNC1(atanh, atanh, 0, + "atanh(x)\n\nReturn the hyperbolic arc tangent (measured in radians) of x.") +FUNC1(ceil, ceil, 0, "ceil(x)\n\nReturn the ceiling of x as a float.\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.") -FUNC1(floor, floor, +FUNC1(floor, floor, 0, "floor(x)\n\nReturn the floor of x as a float.\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 * @@ -172,13 +289,17 @@ 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, @@ -191,19 +312,24 @@ 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, @@ -216,12 +342,10 @@ 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, @@ -260,7 +384,7 @@ } /* Else let libm handle it by itself. */ - return math_1(arg, func); + return math_1(arg, func, 0); } static PyObject * @@ -303,6 +427,141 @@ 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; @@ -356,16 +615,16 @@ "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}, @@ -379,6 +638,7 @@ {"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}, @@ -400,27 +660,15 @@ 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); + PyModule_AddObject(m, "pi", PyFloat_FromDouble(Py_MATH_PI)); + PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E)); - if (!(v = PyFloat_FromDouble(Py_MATH_E))) - goto finally; - if (PyDict_SetItemString(d, "e", v) < 0) - goto finally; - Py_DECREF(v); - - finally: + finally: return; }