From c8e13154c224753f6a21fef1c290cb2f6d3a11f4 Mon Sep 17 00:00:00 2001 From: Pauli Virtanen Date: Thu, 14 Oct 2010 11:04:51 +0200 Subject: [PATCH 1/4] BUG: DOC: fix invalid vdot documentation (cherry picked from commit eca8f94e294003336d901b7e718375fad0c2619c) --- numpy/add_newdocs.py | 23 +++++++++-------------- 1 file changed, 9 insertions(+), 14 deletions(-) diff --git a/numpy/add_newdocs.py b/numpy/add_newdocs.py index cf57c031bfe..3557c00a376 100644 --- a/numpy/add_newdocs.py +++ b/numpy/add_newdocs.py @@ -1295,12 +1295,9 @@ dot(`a`, `b`). If the first argument is complex the complex conjugate of the first argument is used for the calculation of the dot product. - For 2-D arrays it is equivalent to matrix multiplication, and for 1-D - arrays to inner product of vectors (with complex conjugation of `a`). - For N dimensions it is a sum product over the last axis of `a` and - the second-to-last of `b`:: - - dot(a, b)[i,j,k,m] = sum(a[i,j,:] * b[k,:,m]) + Note that `vdot` handles multidimensional arrays differently than `dot`: + it does *not* perform a matrix product, but flattens input arguments + to 1-D vectors first. Consequently, it should only be used for vectors. Parameters ---------- @@ -1313,7 +1310,7 @@ Returns ------- output : ndarray - Returns dot product of `a` and `b`. Can be an int, float, or + Dot product of `a` and `b`. Can be an int, float, or complex depending on the types of `a` and `b`. See Also @@ -1321,13 +1318,6 @@ dot : Return the dot product without using the complex conjugate of the first argument. - Notes - ----- - The dot product is the summation of element wise multiplication. - - .. math:: - a \\cdot b = \\sum_{i=1}^n a_i^*b_i = a_1^*b_1+a_2^*b_2+\\cdots+a_n^*b_n - Examples -------- >>> a = np.array([1+2j,3+4j]) @@ -1336,12 +1326,17 @@ (70-8j) >>> np.vdot(b, a) (70+8j) + + Note that higher-dimensional arrays are flattened! + >>> a = np.array([[1, 4], [5, 6]]) >>> b = np.array([[4, 1], [2, 2]]) >>> np.vdot(a, b) 30 >>> np.vdot(b, a) 30 + >>> 1*4 + 4*1 + 5*2 + 6*2 + 30 """) From 4e177d36755ea5ea70f6512d12f7028639612d3e Mon Sep 17 00:00:00 2001 From: Matthew Brett Date: Mon, 11 Oct 2010 16:09:53 -0700 Subject: [PATCH 2/4] BUG: get fortran arch flags from C arch flags if available. Closes #1399. --- numpy/distutils/fcompiler/gnu.py | 36 ++++++++++---------------------- 1 file changed, 11 insertions(+), 25 deletions(-) diff --git a/numpy/distutils/fcompiler/gnu.py b/numpy/distutils/fcompiler/gnu.py index becdefad57f..632c6d97ad1 100644 --- a/numpy/distutils/fcompiler/gnu.py +++ b/numpy/distutils/fcompiler/gnu.py @@ -202,8 +202,18 @@ def get_flags_opt(self): opt.append('-funroll-loops') return opt + def _c_arch_flags(self): + """ Return detected arch flags from CFLAGS """ + from distutils import sysconfig + cflags = sysconfig.get_config_vars()['CFLAGS'] + arch_re = re.compile(r"-arch\s+(\w+)") + arch_flags = [] + for arch in arch_re.findall(cflags): + arch_flags += ['-arch', arch] + return arch_flags + def get_flags_arch(self): - return [] + return self._c_arch_flags() class Gnu95FCompiler(GnuFCompiler): compiler_type = 'gnu95' @@ -249,30 +259,6 @@ def version_match(self, version_string): g2c = 'gfortran' - def _universal_flags(self, cmd): - """Return a list of -arch flags for every supported architecture.""" - if not sys.platform == 'darwin': - return [] - arch_flags = [] - for arch in ["ppc", "i686", "x86_64", "ppc64"]: - if _can_target(cmd, arch): - arch_flags.extend(["-arch", arch]) - return arch_flags - - def get_flags(self): - flags = GnuFCompiler.get_flags(self) - arch_flags = self._universal_flags(self.compiler_f90) - if arch_flags: - flags[:0] = arch_flags - return flags - - def get_flags_linker_so(self): - flags = GnuFCompiler.get_flags_linker_so(self) - arch_flags = self._universal_flags(self.linker_so) - if arch_flags: - flags[:0] = arch_flags - return flags - def get_library_dirs(self): opt = GnuFCompiler.get_library_dirs(self) if sys.platform == 'win32': From 0e792e60cd32bba6bbcdfaa23af3c379c1ba0dff Mon Sep 17 00:00:00 2001 From: Pauli Virtanen Date: Sat, 16 Oct 2010 11:36:44 +0200 Subject: [PATCH 3/4] BUG: core: adjust ComplexWarning location frame up by one, so that the warning comes from the correct location (#1639) (cherry picked from commit d24db3430a6710d18f9f9e77c89ee442e9c0d631) --- numpy/core/src/multiarray/convert_datatype.c | 2 +- numpy/core/src/scalarmathmodule.c.src | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/numpy/core/src/multiarray/convert_datatype.c b/numpy/core/src/multiarray/convert_datatype.c index d0dd6d271ab..2c66a584b70 100644 --- a/numpy/core/src/multiarray/convert_datatype.c +++ b/numpy/core/src/multiarray/convert_datatype.c @@ -121,7 +121,7 @@ PyArray_GetCastFunc(PyArray_Descr *descr, int type_num) #if PY_VERSION_HEX >= 0x02050000 ret = PyErr_WarnEx(cls, "Casting complex values to real discards the imaginary " - "part", 0); + "part", 1); #else ret = PyErr_Warn(cls, "Casting complex values to real discards the imaginary " diff --git a/numpy/core/src/scalarmathmodule.c.src b/numpy/core/src/scalarmathmodule.c.src index 6261c4c73c4..6913f517dbd 100644 --- a/numpy/core/src/scalarmathmodule.c.src +++ b/numpy/core/src/scalarmathmodule.c.src @@ -990,7 +990,7 @@ emit_complexwarning() #if PY_VERSION_HEX >= 0x02050000 return PyErr_WarnEx(cls, "Casting complex values to real discards the imaginary " - "part", 0); + "part", 1); #else return PyErr_Warn(cls, "Casting complex values to real discards the imaginary " From b55eacd8bf62e592ec44fef351a81c60535d2965 Mon Sep 17 00:00:00 2001 From: Pauli Virtanen Date: Mon, 11 Oct 2010 23:58:50 +0200 Subject: [PATCH 4/4] BUG: core: implement a long-int loop for ldexp, for cases where int != long (#1633) long != int on many 64-bit platforms, so a second ufunc loop is needed to handle ldexp long int inputs. (cherry picked from commit 93f7521dd0ac9edc0034eec5501a126cc4683b70) --- numpy/core/src/umath/loops.c.src | 30 ++++++++++++++++++++++++++ numpy/core/src/umath/loops.h | 6 ++++++ numpy/core/src/umath/loops.h.src | 2 ++ numpy/core/src/umath/umathmodule.c.src | 27 +++++++++++++++++------ numpy/core/tests/test_umath.py | 27 +++++++++++++++++------ 5 files changed, 80 insertions(+), 12 deletions(-) diff --git a/numpy/core/src/umath/loops.c.src b/numpy/core/src/umath/loops.c.src index 24a70cc53ec..df1b5e92104 100644 --- a/numpy/core/src/umath/loops.c.src +++ b/numpy/core/src/umath/loops.c.src @@ -1123,6 +1123,36 @@ NPY_NO_EXPORT void *((@type@ *)op1) = ldexp@c@(in1, in2); } } + +NPY_NO_EXPORT void +@TYPE@_ldexp_long(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func)) +{ + /* + * Additional loop to handle long integer inputs (cf. #866, #1633). + * long != int on many 64-bit platforms, so we need this second loop + * to handle the default integer type. + */ + BINARY_LOOP { + const @type@ in1 = *(@type@ *)ip1; + const long in2 = *(long *)ip2; + if (((int)in2) == in2) { + /* Range OK */ + *((@type@ *)op1) = ldexp@c@(in1, ((int)in2)); + } + else { + /* + * Outside int range -- also ldexp will overflow in this case, + * given that exponent has less bits than int. + */ + if (in2 > 0) { + *((@type@ *)op1) = ldexp@c@(in1, NPY_MAX_INT); + } + else { + *((@type@ *)op1) = ldexp@c@(in1, NPY_MIN_INT); + } + } + } +} #endif #define @TYPE@_true_divide @TYPE@_divide diff --git a/numpy/core/src/umath/loops.h b/numpy/core/src/umath/loops.h index 216179b2537..4234c48d6b7 100644 --- a/numpy/core/src/umath/loops.h +++ b/numpy/core/src/umath/loops.h @@ -1675,6 +1675,8 @@ FLOAT_frexp(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func)); #ifdef HAVE_LDEXPF NPY_NO_EXPORT void FLOAT_ldexp(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func)); +NPY_NO_EXPORT void +FLOAT_ldexp_long(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func)); #endif #define FLOAT_true_divide FLOAT_divide @@ -1827,6 +1829,8 @@ DOUBLE_frexp(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func)) #ifdef HAVE_LDEXP NPY_NO_EXPORT void DOUBLE_ldexp(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func)); +NPY_NO_EXPORT void +DOUBLE_ldexp_long(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func)); #endif #define DOUBLE_true_divide DOUBLE_divide @@ -1979,6 +1983,8 @@ LONGDOUBLE_frexp(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(fu #ifdef HAVE_LDEXPL NPY_NO_EXPORT void LONGDOUBLE_ldexp(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func)); +NPY_NO_EXPORT void +LONGDOUBLE_ldexp_long(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func)); #endif #define LONGDOUBLE_true_divide LONGDOUBLE_divide diff --git a/numpy/core/src/umath/loops.h.src b/numpy/core/src/umath/loops.h.src index 3c277ff4776..4019c388cfc 100644 --- a/numpy/core/src/umath/loops.h.src +++ b/numpy/core/src/umath/loops.h.src @@ -268,6 +268,8 @@ NPY_NO_EXPORT void #ifdef HAVE_LDEXP@C@ NPY_NO_EXPORT void @TYPE@_ldexp(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func)); +NPY_NO_EXPORT void +@TYPE@_ldexp_long(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func)); #endif #define @TYPE@_true_divide @TYPE@_divide diff --git a/numpy/core/src/umath/umathmodule.c.src b/numpy/core/src/umath/umathmodule.c.src index 9694f521d68..c3da9f3d3e0 100644 --- a/numpy/core/src/umath/umathmodule.c.src +++ b/numpy/core/src/umath/umathmodule.c.src @@ -164,6 +164,8 @@ static PyUFuncGenericFunction frexp_functions[] = { }; static void * blank3_data[] = { (void *)NULL, (void *)NULL, (void *)NULL}; +static void * blank6_data[] = { (void *)NULL, (void *)NULL, (void *)NULL, + (void *)NULL, (void *)NULL, (void *)NULL}; static char frexp_signatures[] = { #ifdef HAVE_FREXPF PyArray_FLOAT, PyArray_FLOAT, PyArray_INT, @@ -174,22 +176,35 @@ static char frexp_signatures[] = { #endif }; +#if NPY_SIZEOF_LONG == NPY_SIZEOF_INT +#define LDEXP_LONG(typ) typ##_ldexp +#else +#define LDEXP_LONG(typ) typ##_ldexp_long +#endif + static PyUFuncGenericFunction ldexp_functions[] = { #ifdef HAVE_LDEXPF FLOAT_ldexp, + LDEXP_LONG(FLOAT), #endif - DOUBLE_ldexp + DOUBLE_ldexp, + LDEXP_LONG(DOUBLE) #ifdef HAVE_LDEXPL - ,LONGDOUBLE_ldexp + , + LONGDOUBLE_ldexp, + LDEXP_LONG(LONGDOUBLE) #endif }; static char ldexp_signatures[] = { #ifdef HAVE_LDEXPF PyArray_FLOAT, PyArray_INT, PyArray_FLOAT, + PyArray_FLOAT, PyArray_LONG, PyArray_FLOAT, #endif + PyArray_DOUBLE, PyArray_INT, PyArray_DOUBLE, PyArray_DOUBLE, PyArray_LONG, PyArray_DOUBLE #ifdef HAVE_LDEXPL + ,PyArray_LONGDOUBLE, PyArray_INT, PyArray_LONGDOUBLE ,PyArray_LONGDOUBLE, PyArray_LONG, PyArray_LONGDOUBLE #endif }; @@ -213,14 +228,14 @@ InitOtherOperators(PyObject *dictionary) { PyDict_SetItemString(dictionary, "frexp", f); Py_DECREF(f); - num = 1; + num = 2; #ifdef HAVE_LDEXPL - num += 1; + num += 2; #endif #ifdef HAVE_LDEXPF - num += 1; + num += 2; #endif - f = PyUFunc_FromFuncAndData(ldexp_functions, blank3_data, ldexp_signatures, num, + f = PyUFunc_FromFuncAndData(ldexp_functions, blank6_data, ldexp_signatures, num, 2, 1, PyUFunc_None, "ldexp", "Compute y = x1 * 2**x2.",0); PyDict_SetItemString(dictionary, "ldexp", f); diff --git a/numpy/core/tests/test_umath.py b/numpy/core/tests/test_umath.py index 01b49e5508a..f65a7cd8a82 100644 --- a/numpy/core/tests/test_umath.py +++ b/numpy/core/tests/test_umath.py @@ -387,14 +387,29 @@ def test_nan_any(self): class TestLdexp(TestCase): + def _check_ldexp(self, tp): + assert_almost_equal(ncu.ldexp(np.array(2., np.float32), + np.array(3, tp)), 16.) + assert_almost_equal(ncu.ldexp(np.array(2., np.float64), + np.array(3, tp)), 16.) + assert_almost_equal(ncu.ldexp(np.array(2., np.longdouble), + np.array(3, tp)), 16.) + def test_ldexp(self): + # The default Python int type should work assert_almost_equal(ncu.ldexp(2., 3), 16.) - assert_almost_equal(ncu.ldexp(np.array(2., np.float32), np.array(3, np.int16)), 16.) - assert_almost_equal(ncu.ldexp(np.array(2., np.float32), np.array(3, np.int32)), 16.) - assert_almost_equal(ncu.ldexp(np.array(2., np.float64), np.array(3, np.int16)), 16.) - assert_almost_equal(ncu.ldexp(np.array(2., np.float64), np.array(3, np.int32)), 16.) - assert_almost_equal(ncu.ldexp(np.array(2., np.longdouble), np.array(3, np.int16)), 16.) - assert_almost_equal(ncu.ldexp(np.array(2., np.longdouble), np.array(3, np.int32)), 16.) + # The following int types should all be accepted + self._check_ldexp(np.int8) + self._check_ldexp(np.int16) + self._check_ldexp(np.int32) + self._check_ldexp('i') + self._check_ldexp('l') + + def test_ldexp_overflow(self): + imax = np.iinfo(np.dtype('l')).max + imin = np.iinfo(np.dtype('l')).min + assert_equal(ncu.ldexp(2., imax), np.inf) + assert_equal(ncu.ldexp(2., imin), 0) class TestMaximum(TestCase):