diff options
70 files changed, 3020 insertions, 1061 deletions
diff --git a/.gitignore b/.gitignore index 8387f4a03..ceae918fc 100644 --- a/.gitignore +++ b/.gitignore @@ -6,6 +6,10 @@ *$ *.bak *.diff +.idea/ +*.iml +*.ipr +*.iws *.org .project *.rej diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md new file mode 100644 index 000000000..937fbe090 --- /dev/null +++ b/CONTRIBUTING.md @@ -0,0 +1,26 @@ +# Contributing to numpy + +## Reporting issues + +When reporting issues please include as much detail as possible about your +operating system, numpy version and python version. Whenever possible, please +also include a brief, self-contained code example that demonstrates the problem. + +If you are reporting a segfault please include a GDB traceback, which you can +generate by following +[these instructions.](https://github.com/numpy/numpy/blob/master/doc/source/dev/development_environment.rst#debugging) + +## Contributing code + +Thanks for your interest in contributing code to numpy! + ++ If this is your first time contributing to a project on GitHub, please read +through our +[guide to contributing to numpy](http://docs.scipy.org/doc/numpy-dev/dev/index.html) ++ If you have contributed to other projects on GitHub you can go straight to our +[development workflow](http://docs.scipy.org/doc/numpy-dev/dev/gitwash/development_workflow.html) + +Either way, please be sure to follow our +[convention for commit messages](http://docs.scipy.org/doc/numpy-dev/dev/gitwash/development_workflow.html). + + diff --git a/INSTALL.txt b/INSTALL.txt index 9c12ba602..01726df04 100644 --- a/INSTALL.txt +++ b/INSTALL.txt @@ -28,7 +28,7 @@ Building NumPy requires the following software installed: Python must also be compiled with the zlib module enabled. -2) nose__ (optional) 0.10.3 or later +2) nose__ (optional) 1.0 or later This is required for testing numpy, but not for using it. diff --git a/bento.info b/bento.info index 2348b44be..29bde06cc 100644 --- a/bento.info +++ b/bento.info @@ -38,8 +38,6 @@ ExtraSourceFiles: setup.py, numpy/core/include/numpy/*.h, numpy/core/include/numpy/*.in, - numpy/core/include/numpy/fenv/*.h, - numpy/core/include/numpy/fenv/*.c, numpy/core/*.ini.in, numpy/core/src/npymath/*.h, numpy/core/src/multiarray/*.c, @@ -73,8 +71,7 @@ DataFiles: f2py-data DataFiles: numpy-includes TargetDir: $sitedir Files: - numpy/core/include/numpy/*.h, - numpy/core/include/numpy/fenv/*.h + numpy/core/include/numpy/*.h HookFile: bscript MetaTemplateFiles: numpy/version.py.in, numpy/__config__.py.in diff --git a/doc/HOWTO_DOCUMENT.rst.txt b/doc/HOWTO_DOCUMENT.rst.txt index eac97512a..de0190084 100644 --- a/doc/HOWTO_DOCUMENT.rst.txt +++ b/doc/HOWTO_DOCUMENT.rst.txt @@ -13,11 +13,15 @@ A Guide to NumPy/SciPy Documentation numpy conventions, you should use the ``numpydoc`` extension so that your docstrings will be handled correctly. For example, Sphinx will extract the ``Parameters`` section from your docstring and convert it into a field - list. Using numpydoc will also avoid the reStructuredText errors produced + list. Using ``numpydoc`` will also avoid the reStructuredText errors produced by plain Sphinx when it encounters numpy docstring conventions like section headers (e.g. ``-------------``) that sphinx does not expect to find in docstrings. + Some features described in this document require a recent version of + ``numpydoc``. For example, the **Yields** section was added in + ``numpydoc`` 0.6. + It is available from: * `numpydoc on PyPI <http://pypi.python.org/pypi/numpydoc>`_ @@ -163,14 +167,14 @@ The sections of the docstring are: `ndobj_old` will be removed in Numpy 2.0, it is replaced by `ndobj_new` because the latter works also with array subclasses. -3. **Extended summary** +3. **Extended Summary** A few sentences giving an extended description. This section should be used to clarify *functionality*, not to discuss implementation detail or background theory, which should rather be - explored in the **notes** section below. You may refer to the + explored in the **Notes** section below. You may refer to the parameters and the function name, but parameter descriptions still - belong in the **parameters** section. + belong in the **Parameters** section. 4. **Parameters** @@ -228,7 +232,7 @@ The sections of the docstring are: 5. **Returns** Explanation of the returned values and their types. Similar to the - **parameters** section, except the name of each return value is optional. + **Parameters** section, except the name of each return value is optional. The type of each return value is always required:: Returns @@ -236,8 +240,8 @@ The sections of the docstring are: int Description of anonymous integer return value. - If both the name and type are specified, the **returns** section takes the - same form as the **parameters** section:: + If both the name and type are specified, the **Returns** section takes the + same form as the **Parameters** section:: Returns ------- @@ -246,13 +250,37 @@ The sections of the docstring are: err_msg : str or None Human readable error message, or None on success. -6. **Other parameters** +6. **Yields** + + Explanation of the yielded values and their types. This is relevant to + generators only. Similar to the **Returns** section in that the name of + each value is optional, but the type of each value is always required:: + + Yields + ------ + int + Description of the anonymous integer return value. + + If both the name and type are specified, the **Yields** section takes the + same form as the **Returns** section:: + + Yields + ------ + err_code : int + Non-zero value indicates error code, or zero on success. + err_msg : str or None + Human readable error message, or None on success. + + Support for the **Yields** section was added in `numpydoc + <https://github.com/numpy/numpydoc>`_ version 0.6. + +7. **Other Parameters** An optional section used to describe infrequently used parameters. It should only be used if a function has a large number of keyword - parameters, to prevent cluttering the **parameters** section. + parameters, to prevent cluttering the **Parameters** section. -7. **Raises** +8. **Raises** An optional section detailing which errors get raised and under what conditions:: @@ -265,7 +293,7 @@ The sections of the docstring are: This section should be used judiciously, i.e., only for errors that are non-obvious or have a large chance of getting raised. -8. **See Also** +9. **See Also** An optional section used to refer to related code. This section can be very useful, but should be used judiciously. The goal is to @@ -304,7 +332,7 @@ The sections of the docstring are: func_b, func_c_, func_d func_e -9. **Notes** +10. **Notes** An optional section that provides additional information about the code, possibly including a discussion of the algorithm. This @@ -349,7 +377,7 @@ The sections of the docstring are: where filename is a path relative to the reference guide source directory. -10. **References** +11. **References** References cited in the **notes** section may be listed here, e.g. if you cited the article below using the text ``[1]_``, @@ -374,7 +402,7 @@ The sections of the docstring are: should not be required to understand it. References are numbered, starting from one, in the order in which they are cited. -11. **Examples** +12. **Examples** An optional section for examples, using the `doctest <http://docs.python.org/library/doctest.html>`_ format. @@ -437,10 +465,10 @@ Class docstring ``````````````` Use the same sections as outlined above (all except ``Returns`` are applicable). The constructor (``__init__``) should also be documented -here, the **parameters** section of the docstring details the constructors +here, the **Parameters** section of the docstring details the constructors parameters. -An **Attributes** section, located below the **parameters** section, +An **Attributes** section, located below the **Parameters** section, may be used to describe class variables:: Attributes @@ -466,7 +494,7 @@ In general, it is not necessary to list class methods. Those that are not part of the public API have names that start with an underscore. In some cases, however, a class may have a great many methods, of which only a few are relevant (e.g., subclasses of ndarray). Then, it -becomes useful to have an additional **methods** section:: +becomes useful to have an additional **Methods** section:: class Photo(ndarray): """ @@ -489,8 +517,8 @@ becomes useful to have an additional **methods** section:: """ If it is necessary to explain a private method (use with care!), it can -be referred to in the **extended summary** or the **notes**. Do not -list private methods in the **methods** section. +be referred to in the **Extended Summary** or the **Notes** section. +Do not list private methods in the **methods** section. Note that `self` is *not* listed as the first parameter of methods. @@ -501,7 +529,8 @@ Document these as you would any other function. Do not include (which is the case for many ndarray methods for example), the function docstring should contain the detailed documentation, and the method docstring should refer to it. Only put brief summary and **See Also** sections in the -method docstring. +method docstring. The method should use a **Returns** or **Yields** section, +as appropriate. Documenting class instances @@ -520,6 +549,14 @@ instances a useful docstring, we do the following: sections. +Documenting generators +---------------------- +Generators should be documented just as functions are documented. The +only difference is that one should use the **Yields** section instead +of the **Returns** section. Support for the **Yields** section was added in +`numpydoc <https://github.com/numpy/numpydoc>`_ version 0.6. + + Documenting constants --------------------- Use the same sections as outlined for functions where applicable:: diff --git a/doc/neps/ufunc-overrides.rst b/doc/neps/ufunc-overrides.rst index 5a0a0334f..be3b8b4aa 100644 --- a/doc/neps/ufunc-overrides.rst +++ b/doc/neps/ufunc-overrides.rst @@ -3,7 +3,7 @@ A Mechanism for Overriding Ufuncs ================================= :Author: Blake Griffith -:Contact: blake.g@utexa.edu +:Contact: blake.g@utexas.edu :Date: 2013-07-10 :Author: Pauli Virtanen diff --git a/doc/release/1.10.0-notes.rst b/doc/release/1.10.0-notes.rst index 76538760e..f9d202ec3 100644 --- a/doc/release/1.10.0-notes.rst +++ b/doc/release/1.10.0-notes.rst @@ -61,6 +61,15 @@ C API The changes to *swapaxes* also apply to the *PyArray_SwapAxes* C function, which now returns a view in all cases. +recarray field return types +~~~~~~~~~~~~~~~~~~~~~~~~~~~ +Previously the returned types for recarray fields accessed by attribute and by +index were inconsistent, and fields of string type were returned as chararrays. +Now, fields accessed by either attribute or indexing will return an ndarray for +fields of non-structured type, and a recarray for fields of structured type. +Notably, this affect recarrays containing strings with whitespace, as trailing +whitespace is trimmed from chararrays but kept in ndarrays of string type. +Also, the dtype.type of nested structured fields is now inherited. New Features ============ @@ -78,6 +87,13 @@ extensions is now performed in *n* parallel processes. The parallelization is limited to files within one extension so projects using Cython will not profit because it builds extensions from single files. +*genfromtxt* has an new ``max_rows`` argument +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +A ``max_rows`` argument has been added to *genfromtxt* to limit the +number of rows read in a single call. Using this functionality, it is +possible to read in multiple arrays stored in a single file by making +repeated calls to the function. + Improvements ============ @@ -105,11 +121,41 @@ close the previous and the next period cycles, resulting in the correct interpolation behavior. *np.pad* supports more input types for ``pad_width`` and ``constant_values`` -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ``constant_values`` parameters now accepts NumPy arrays and float values. NumPy arrays are supported as input for ``pad_width``, and an exception is raised if its values are not of integral type. +More system C99 complex functions detected and used +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +All of the functions ``in complex.h`` are now detected. There are new +fallback implementations of the following functions. + +* npy_ctan, +* npy_cacos, npy_casin, npy_catan +* npy_ccosh, npy_csinh, npy_ctanh, +* npy_cacosh, npy_casinh, npy_catanh + +As a result of these improvements, there will be some small changes in +returned values, especially for corner cases. + +*np.loadtxt* support for the strings produced by the ``float.hex`` method +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +The strings produced by ``float.hex`` look like ``0x1.921fb54442d18p+1``, +so this is not the hex used to represent unsigned integer types. + +*np.isclose* properly handles minimal values of integer dtypes +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +In order to properly handle minimal values of integer types, *np.isclose* will +now cast to the float dtype during comparisons. This aligns its behavior with +what was provided by *np.allclose*. + +*np.allclose* uses *np.isclose* internally. +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +*np.allcose* now uses *np.isclose* internally and inherits the ability to +compare NaNs as equal by setting ``equal_nan=True``. Subclasses, such as +*np.ma.MaskedArray*, are also preserved now. + Changes ======= diff --git a/doc/source/dev/development_environment.rst b/doc/source/dev/development_environment.rst new file mode 100644 index 000000000..6221353ce --- /dev/null +++ b/doc/source/dev/development_environment.rst @@ -0,0 +1,211 @@ +.. _development-environment: + +Setting up and using your development environment +================================================= + + +Recommended development setup +----------------------------- + +Since NumPy contains parts written in C and Cython that need to be +compiled before use, make sure you have the necessary compilers and Python +development headers installed - see :ref:`building-from-source`. + +Having compiled code also means that importing Numpy from the development +sources needs some additional steps, which are explained below. For the rest +of this chapter we assume that you have set up your git repo as described in +:ref:`using-git`. + +To build the development version of NumPy and run tests, spawn +interactive shells with the Python import paths properly set up etc., +do one of:: + + $ python runtests.py -v + $ python runtests.py -v -s random + $ python runtests.py -v -t numpy/core/tests/test_iter.py:test_iter_c_order + $ python runtests.py --ipython + $ python runtests.py --python somescript.py + $ python runtests.py --bench + $ python runtests.py -g -m full + +This builds Numpy first, so the first time it may take a few minutes. If +you specify ``-n``, the tests are run against the version of NumPy (if +any) found on current PYTHONPATH. + +Using ``runtests.py`` is the recommended approach to running tests. +There are also a number of alternatives to it, for example in-place +build or installing to a virtualenv. See the FAQ below for details. + + +Building in-place +----------------- + +For development, you can set up an in-place build so that changes made to +``.py`` files have effect without rebuild. First, run:: + + $ python setup.py build_ext -i + +This allows you to import the in-place built NumPy *from the repo base +directory only*. If you want the in-place build to be visible outside that +base dir, you need to point your ``PYTHONPATH`` environment variable to this +directory. Some IDEs (Spyder for example) have utilities to manage +``PYTHONPATH``. On Linux and OSX, you can run the command:: + + $ export PYTHONPATH=$PWD + +and on Windows:: + + $ set PYTHONPATH=/path/to/numpy + +Now editing a Python source file in NumPy allows you to immediately +test and use your changes (in ``.py`` files), by simply restarting the +interpreter. + +Note that another way to do an inplace build visible outside the repo base dir +is with ``python setup.py develop``. This doesn't work for NumPy, because +NumPy builds don't use ``setuptools`` by default. ``python setupegg.py +develop`` will work though. + + +Other build options +------------------- + +It's possible to do a parallel build with ``numpy.distutils`` with the ``-j`` option; +see :ref:`parallel-builds` for more details. + +In order to install the development version of NumPy in ``site-packages``, use +``python setup.py install --user``. + +A similar approach to in-place builds and use of ``PYTHONPATH`` but outside the +source tree is to use:: + + $ python setup.py install --prefix /some/owned/folder + $ export PYTHONPATH=/some/owned/folder/lib/python3.4/site-packages + +Besides ``numpy.distutils``, NumPy supports building with `Bento`_. +This provides (among other things) faster builds and a build log that's much +more readable than the ``distutils`` one. Note that support is still fairly +experimental, partly due to Bento relying on `Waf`_ which tends to have +non-backwards-compatible API changes. Working versions of Bento and Waf are +run on TravisCI, see ``tools/travis-test.sh``. + + +Using virtualenvs +----------------- + +A frequently asked question is "How do I set up a development version of NumPy +in parallel to a released version that I use to do my job/research?". + +One simple way to achieve this is to install the released version in +site-packages, by using a binary installer or pip for example, and set +up the development version in a virtualenv. First install +`virtualenv`_ (optionally use `virtualenvwrapper`_), then create your +virtualenv (named numpy-dev here) with:: + + $ virtualenv numpy-dev + +Now, whenever you want to switch to the virtual environment, you can use the +command ``source numpy-dev/bin/activate``, and ``deactivate`` to exit from the +virtual environment and back to your previous shell. + + +Running tests +------------- + +Besides using ``runtests.py``, there are various ways to run the tests. Inside +the interpreter, tests can be run like this:: + + >>> np.test() + >>> np.test('full') # Also run tests marked as slow + >>> np.test('full', verbose=2) # Additionally print test name/file + +Or a similar way from the command line:: + + $ python -c "import numpy as np; np.test()" + +Tests can also be run with ``nosetests numpy``, however then the NumPy-specific +``nose`` plugin is not found which causes tests marked as ``KnownFailure`` to +be reported as errors. + +Running individual test files can be useful; it's much faster than running the +whole test suite or that of a whole module (example: ``np.random.test()``). +This can be done with:: + + $ python path_to_testfile/test_file.py + +That also takes extra arguments, like ``--pdb`` which drops you into the Python +debugger when a test fails or an exception is raised. + +Running tests with `tox`_ is also supported. For example, to build NumPy and +run the test suite with Python 3.4, use:: + + $ tox -e py34 + +For more extensive info on running and writing tests, see +https://github.com/numpy/numpy/blob/master/doc/TESTS.rst.txt . + + +Rebuilding & cleaning the workspace +----------------------------------- + +Rebuilding NumPy after making changes to compiled code can be done with the +same build command as you used previously - only the changed files will be +re-built. Doing a full build, which sometimes is necessary, requires cleaning +the workspace first. The standard way of doing this is (*note: deletes any +uncommitted files!*):: + + $ git clean -xdf + +When you want to discard all changes and go back to the last commit in the +repo, use one of:: + + $ git checkout . + $ git reset --hard + + +Debugging +--------- + +Another frequently asked question is "How do I debug C code inside Numpy?". +The easiest way to do this is to first write a Python script that invokes the C +code whose execution you want to debug. For instance ``mytest.py``:: + + from numpy import linspace + x = np.arange(5) + np.empty_like(x) + +Now, you can run:: + + $ gdb --args python runtests.py -g --python mytest.py + +And then in the debugger:: + + (gdb) break array_empty_like + (gdb) run + +The execution will now stop at the corresponding C function and you can step +through it as usual. With the Python extensions for gdb installed (often the +default on Linux), a number of useful Python-specific commands are available. +For example to see where in the Python code you are, use ``py-list``. For more +details, see `DebuggingWithGdb`_. + +Instead of plain ``gdb`` you can of course use your favourite +alternative debugger; run it on the python binary with arguments +``runtests.py -g --python mytest.py``. + +Building NumPy with a Python built with debug support (on Linux distributions +typically packaged as ``python-dbg``) is highly recommended. + + + +.. _Bento: http://cournape.github.io/Bento/ + +.. _DebuggingWithGdb: https://wiki.python.org/moin/DebuggingWithGdb + +.. _tox: http://tox.testrun.org + +.. _virtualenv: http://www.virtualenv.org/ + +.. _virtualenvwrapper: http://www.doughellmann.com/projects/virtualenvwrapper/ + +.. _Waf: https://code.google.com/p/waf/ diff --git a/doc/source/dev/gitwash/index.rst b/doc/source/dev/gitwash/index.rst index 9d733dd1c..ae7ce69de 100644 --- a/doc/source/dev/gitwash/index.rst +++ b/doc/source/dev/gitwash/index.rst @@ -1,7 +1,7 @@ .. _using-git: Working with *NumPy* source code -====================================== +================================ Contents: diff --git a/doc/source/dev/index.rst b/doc/source/dev/index.rst index 2229f3ccb..b0d0ec483 100644 --- a/doc/source/dev/index.rst +++ b/doc/source/dev/index.rst @@ -6,5 +6,6 @@ Contributing to Numpy :maxdepth: 3 gitwash/index + development_environment For core developers: see :ref:`development-workflow`. diff --git a/doc/source/reference/c-api.array.rst b/doc/source/reference/c-api.array.rst index 9b8cc04b6..f5f753292 100644 --- a/doc/source/reference/c-api.array.rst +++ b/doc/source/reference/c-api.array.rst @@ -593,18 +593,16 @@ From other objects .. cfunction:: PyObject* PyArray_FromStructInterface(PyObject* op) Returns an ndarray object from a Python object that exposes the - :obj:`__array_struct__`` method and follows the array interface - protocol. If the object does not contain this method then a + :obj:`__array_struct__` attribute and follows the array interface + protocol. If the object does not contain this attribute then a borrowed reference to :cdata:`Py_NotImplemented` is returned. .. cfunction:: PyObject* PyArray_FromInterface(PyObject* op) Returns an ndarray object from a Python object that exposes the - :obj:`__array_shape__` and :obj:`__array_typestr__` - methods following - the array interface protocol. If the object does not contain one - of these method then a borrowed reference to :cdata:`Py_NotImplemented` - is returned. + :obj:`__array_interface__` attribute following the array interface + protocol. If the object does not contain this attribute then a + borrowed reference to :cdata:`Py_NotImplemented` is returned. .. cfunction:: PyObject* PyArray_FromArrayAttr(PyObject* op, PyArray_Descr* dtype, PyObject* context) @@ -806,15 +804,28 @@ General check of Python Type sub-type of :cdata:`PyGenericArr_Type` ), or an instance of (a sub-class of) :cdata:`PyArray_Type` whose dimensionality is 0. +.. cfunction:: PyArray_IsPythonNumber(op) + + Evaluates true if *op* is an instance of a builtin numeric type (int, + float, complex, long, bool) + .. cfunction:: PyArray_IsPythonScalar(op) - Evaluates true if *op* is a builtin Python "scalar" object (int, + Evaluates true if *op* is a builtin Python scalar object (int, float, complex, str, unicode, long, bool). .. cfunction:: PyArray_IsAnyScalar(op) - Evaluates true if *op* is either a Python scalar or an array - scalar (an instance of a sub- type of :cdata:`PyGenericArr_Type` ). + Evaluates true if *op* is either a Python scalar object (see + :cfunc:`PyArray_IsPythonScalar`) or an array scalar (an instance of a sub- + type of :cdata:`PyGenericArr_Type` ). + +.. cfunction:: PyArray_CheckAnyScalar(op) + + Evaluates true if *op* is a Python scalar object (see + :cfunc:`PyArray_IsPythonScalar`), an array scalar (an instance of a + sub-type of :cdata:`PyGenericArr_Type`) or an instance of a sub-type of + :cdata:`PyArray_Type` whose dimensionality is 0. Data-type checking @@ -2519,6 +2530,8 @@ Array Scalars .. cfunction:: PyObject* PyArray_Return(PyArrayObject* arr) + This function steals a reference to *arr*. + This function checks to see if *arr* is a 0-dimensional array and, if so, returns the appropriate array scalar. It should be used whenever 0-dimensional arrays could be returned to Python. diff --git a/doc/source/reference/c-api.types-and-structures.rst b/doc/source/reference/c-api.types-and-structures.rst index 473e25010..43abe24c7 100644 --- a/doc/source/reference/c-api.types-and-structures.rst +++ b/doc/source/reference/c-api.types-and-structures.rst @@ -407,7 +407,10 @@ PyArrayDescr_Type PyArray_ScalarKindFunc *scalarkind; int **cancastscalarkindto; int *cancastto; - int listpickle + PyArray_FastClipFunc *fastclip; + PyArray_FastPutmaskFunc *fastputmask; + PyArray_FastTakeFunc *fasttake; + PyArray_ArgFunc *argmin; } PyArray_ArrFuncs; The concept of a behaved segment is used in the description of the @@ -417,8 +420,7 @@ PyArrayDescr_Type functions can (and must) deal with mis-behaved arrays. The other functions require behaved memory segments. - .. cmember:: void cast(void *from, void *to, npy_intp n, void *fromarr, - void *toarr) + .. cmember:: void cast(void *from, void *to, npy_intp n, void *fromarr, void *toarr) An array of function pointers to cast from the current type to all of the other builtin types. Each function casts a @@ -444,8 +446,7 @@ PyArrayDescr_Type a zero is returned, otherwise, a negative one is returned (and a Python error set). - .. cmember:: void copyswapn(void *dest, npy_intp dstride, void *src, - npy_intp sstride, npy_intp n, int swap, void *arr) + .. cmember:: void copyswapn(void *dest, npy_intp dstride, void *src, npy_intp sstride, npy_intp n, int swap, void *arr) .. cmember:: void copyswap(void *dest, void *src, int swap, void *arr) @@ -471,8 +472,7 @@ PyArrayDescr_Type ``d2``, and -1 if * ``d1`` < * ``d2``. The array object ``arr`` is used to retrieve itemsize and field information for flexible arrays. - .. cmember:: int argmax(void* data, npy_intp n, npy_intp* max_ind, - void* arr) + .. cmember:: int argmax(void* data, npy_intp n, npy_intp* max_ind, void* arr) A pointer to a function that retrieves the index of the largest of ``n`` elements in ``arr`` beginning at the element @@ -481,8 +481,7 @@ PyArrayDescr_Type always 0. The index of the largest element is returned in ``max_ind``. - .. cmember:: void dotfunc(void* ip1, npy_intp is1, void* ip2, npy_intp is2, - void* op, npy_intp n, void* arr) + .. cmember:: void dotfunc(void* ip1, npy_intp is1, void* ip2, npy_intp is2, void* op, npy_intp n, void* arr) A pointer to a function that multiplies two ``n`` -length sequences together, adds them, and places the result in @@ -532,8 +531,7 @@ PyArrayDescr_Type computed by repeatedly adding this computed delta. The data buffer must be well-behaved. - .. cmember:: void fillwithscalar(void* buffer, npy_intp length, - void* value, void* arr) + .. cmember:: void fillwithscalar(void* buffer, npy_intp length, void* value, void* arr) A pointer to a function that fills a contiguous ``buffer`` of the given ``length`` with a single scalar ``value`` whose @@ -548,14 +546,13 @@ PyArrayDescr_Type :cdata:`NPY_MERGESORT` are defined). These sorts are done in-place assuming contiguous and aligned data. - .. cmember:: int argsort(void* start, npy_intp* result, npy_intp length, - void \*arr) + .. cmember:: int argsort(void* start, npy_intp* result, npy_intp length, void *arr) An array of function pointers to sorting algorithms for this data type. The same sorting algorithms as for sort are available. The indices producing the sort are returned in - result (which must be initialized with indices 0 to length-1 - inclusive). + ``result`` (which must be initialized with indices 0 to + ``length-1`` inclusive). .. cmember:: PyObject *castdict @@ -587,9 +584,48 @@ PyArrayDescr_Type can be cast to safely (this usually means without losing precision). - .. cmember:: int listpickle + .. cmember:: void fastclip(void *in, npy_intp n_in, void *min, void *max, void *out) + + A function that reads ``n_in`` items from ``in``, and writes to + ``out`` the read value if it is within the limits pointed to by + ``min`` and ``max``, or the corresponding limit if outside. The + memory segments must be contiguous and behaved, and either + ``min`` or ``max`` may be ``NULL``, but not both. + + .. cmember:: void fastputmask(void *in, void *mask, npy_intp n_in, void *values, npy_intp nv) + + A function that takes a pointer ``in`` to an array of ``n_in`` + items, a pointer ``mask`` to an array of ``n_in`` boolean + values, and a pointer ``vals`` to an array of ``nv`` items. + Items from ``vals`` are copied into ``in`` wherever the value + in ``mask`` is non-zero, tiling ``vals`` as needed if + ``nv < n_in``. All arrays must be contiguous and behaved. + + .. cmember:: void fasttake(void *dest, void *src, npy_intp *indarray, npy_intp nindarray, npy_intp n_outer, npy_intp m_middle, npy_intp nelem, NPY_CLIPMODE clipmode) + + A function that takes a pointer ``src`` to a C contiguous, + behaved segment, interpreted as a 3-dimensional array of shape + ``(n_outer, nindarray, nelem)``, a pointer ``indarray`` to a + contiguous, behaved segment of ``m_middle`` integer indices, + and a pointer ``dest`` to a C contiguous, behaved segment, + interpreted as a 3-dimensional array of shape + ``(n_outer, m_middle, nelem)``. The indices in ``indarray`` are + used to index ``src`` along the second dimension, and copy the + corresponding chunks of ``nelem`` items into ``dest``. + ``clipmode`` (which can take on the values :cdata:`NPY_RAISE`, + :cdata:`NPY_WRAP` or :cdata:`NPY_CLIP`) determines how will + indices smaller than 0 or larger than ``nindarray`` will be + handled. + + .. cmember:: int argmin(void* data, npy_intp n, npy_intp* min_ind, void* arr) + + A pointer to a function that retrieves the index of the + smallest of ``n`` elements in ``arr`` beginning at the element + pointed to by ``data``. This function requires that the + memory segment be contiguous and behaved. The return value is + always 0. The index of the smallest element is returned in + ``min_ind``. - Unused. The :cdata:`PyArray_Type` typeobject implements many of the features of Python objects including the tp_as_number, tp_as_sequence, diff --git a/doc/source/reference/swig.testing.rst b/doc/source/reference/swig.testing.rst index decc681c5..c0daaec66 100644 --- a/doc/source/reference/swig.testing.rst +++ b/doc/source/reference/swig.testing.rst @@ -10,8 +10,8 @@ data types are supported, each with 74 different argument signatures, for a total of 888 typemaps supported "out of the box". Each of these typemaps, in turn, might require several unit tests in order to verify expected behavior for both proper and improper inputs. Currently, -this results in 1,427 individual unit tests that are performed when -``make test`` is run in the ``numpy/docs/swig`` subdirectory. +this results in more than 1,000 individual unit tests executed when +``make test`` is run in the ``numpy/tools/swig`` subdirectory. To facilitate this many similar unit tests, some high-level programming techniques are employed, including C and `SWIG`_ macros, diff --git a/doc/source/user/basics.io.genfromtxt.rst b/doc/source/user/basics.io.genfromtxt.rst index edf48bc15..11205e555 100644 --- a/doc/source/user/basics.io.genfromtxt.rst +++ b/doc/source/user/basics.io.genfromtxt.rst @@ -94,12 +94,12 @@ This behavior can be overwritten by setting the optional argument >>> data = "1, abc , 2\n 3, xxx, 4" >>> # Without autostrip - >>> np.genfromtxt(StringIO(data), dtype="|S5") + >>> np.genfromtxt(StringIO(data), delimiter=",", dtype="|S5") array([['1', ' abc ', ' 2'], ['3', ' xxx', ' 4']], dtype='|S5') >>> # With autostrip - >>> np.genfromtxt(StringIO(data), dtype="|S5", autostrip=True) + >>> np.genfromtxt(StringIO(data), delimiter=",", dtype="|S5", autostrip=True) array([['1', 'abc', '2'], ['3', 'xxx', '4']], dtype='|S5') diff --git a/doc/source/user/install.rst b/doc/source/user/install.rst index 29aeff6a3..dcf20498c 100644 --- a/doc/source/user/install.rst +++ b/doc/source/user/install.rst @@ -12,15 +12,15 @@ Windows ------- Good solutions for Windows are, `Enthought Canopy -<https://www.enthought.com/products/canopy/>`_ (which provides binary -installers for Windows, OS X and Linux) and `Python (x, y) -<http://www.pythonxy.com>`_. Both of these packages include Python, NumPy and -many additional packages. +<https://www.enthought.com/products/canopy/>`_, `Anaconda +<http://continuum.io/downloads.html>`_ (which both provide binary installers +for Windows, OS X and Linux) and `Python (x, y) <http://www.pythonxy.com>`_. +Both of these packages include Python, NumPy and many additional packages. A lightweight alternative is to download the Python installer from `www.python.org <http://www.python.org>`_ and the NumPy -installer for your Python version from the Sourceforge `download site <http:// -sourceforge.net/project/showfiles.php?group_id=1369&package_id=175103>`_ +installer for your Python version from the Sourceforge `download site +<http://sourceforge.net/projects/numpy/files/NumPy`_. The NumPy installer includes binaries for different CPU's (without SSE instructions, with SSE2 or with SSE3) and installs the correct one @@ -28,25 +28,27 @@ automatically. If needed, this can be bypassed from the command line with :: numpy-<1.y.z>-superpack-win32.exe /arch nosse -or 'sse2' or 'sse3' instead of 'nosse'. +or ``sse2`` or ``sse3`` instead of ``nosse``. Linux ----- -Most of the major distributions provide packages for NumPy, but these can lag -behind the most recent NumPy release. Pre-built binary packages for Ubuntu are -available on the `scipy ppa -<https://edge.launchpad.net/~scipy/+archive/ppa>`_. Redhat binaries are -available in the `Enthought Canopy -<https://www.enthought.com/products/canopy/>`_. +All major distributions provide packages for NumPy. These are usually +reasonably up-to-date, but sometimes lag behind the most recent NumPy release. Mac OS X -------- -A universal binary installer for NumPy is available from the `download site -<http://sourceforge.net/project/showfiles.php?group_id=1369& -package_id=175103>`_. The `Enthought Canopy -<https://www.enthought.com/products/canopy/>`_ provides NumPy binaries. +Universal binary installers for NumPy are available from the `download site +<http://sourceforge.net/projects/numpy/files/NumPy`_, and wheel packages +from PyPi. With a recent version of `pip`<https://pip.pypa.io/en/latest/>`_ +this will give you a binary install (from the wheel packages) compatible with +at python.org Python, Homebrew and MacPorts:: + + pip install numpy + + +.. _building-from-source: Building from source ==================== @@ -82,7 +84,7 @@ Building NumPy requires the following software installed: Note that NumPy is developed mainly using GNU compilers. Compilers from other vendors such as Intel, Absoft, Sun, NAG, Compaq, Vast, Porland, Lahey, HP, IBM, Microsoft are only supported in the form of community - feedback, and may not work out of the box. GCC 3.x (and later) compilers + feedback, and may not work out of the box. GCC 4.x (and later) compilers are recommended. 3) Linear Algebra libraries @@ -93,6 +95,42 @@ Building NumPy requires the following software installed: can be used, including optimized LAPACK libraries such as ATLAS, MKL or the Accelerate/vecLib framework on OS X. +Basic Installation +------------------ + +To install NumPy run:: + + python setup.py install + +To perform an in-place build that can be run from the source folder run:: + + python setup.py build_ext --inplace + +The NumPy build system uses ``distutils`` and ``numpy.distutils``. +``setuptools`` is only used when building via ``pip`` or with ``python +setupegg.py``. Using ``virtualenv`` should work as expected. + +*Note: for build instructions to do development work on NumPy itself, see +:ref:`development-environment`*. + +.. _parallel-builds: + +Parallel builds +~~~~~~~~~~~~~~~ + +From NumPy 1.10.0 on it's also possible to do a parallel build with:: + + python setup.py build -j 4 install --prefix $HOME/.local + +This will compile numpy on 4 CPUs and install it into the specified prefix. +to perform a parallel in-place build, run:: + + python setup.py build_ext --inplace -j 4 + +The number of build jobs can also be specified via the environment variable +``NPY_NUM_BUILD_JOBS``. + + FORTRAN ABI mismatch -------------------- @@ -147,35 +185,10 @@ Additional compiler flags can be supplied by setting the ``OPT``, Building with ATLAS support --------------------------- -Ubuntu 8.10 (Intrepid) and 9.04 (Jaunty) -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +Ubuntu +~~~~~~ -You can install the necessary packages for optimized ATLAS with this command:: +You can install the necessary package for optimized ATLAS with this command:: sudo apt-get install libatlas-base-dev -If you have a recent CPU with SIMD suppport (SSE, SSE2, etc...), you should -also install the corresponding package for optimal performances. For example, -for SSE2:: - - sudo apt-get install libatlas3gf-sse2 - -This package is not available on amd64 platforms. - -*NOTE*: Ubuntu changed its default fortran compiler from g77 in Hardy to -gfortran in Intrepid. If you are building ATLAS from source and are upgrading -from Hardy to Intrepid or later versions, you should rebuild everything from -scratch, including lapack. - -Ubuntu 8.04 and lower -~~~~~~~~~~~~~~~~~~~~~ - -You can install the necessary packages for optimized ATLAS with this command:: - - sudo apt-get install atlas3-base-dev - -If you have a recent CPU with SIMD suppport (SSE, SSE2, etc...), you should -also install the corresponding package for optimal performances. For example, -for SSE2:: - - sudo apt-get install atlas3-sse2 diff --git a/numpy/add_newdocs.py b/numpy/add_newdocs.py index 72aaf5ec7..7dd8c5649 100644 --- a/numpy/add_newdocs.py +++ b/numpy/add_newdocs.py @@ -664,13 +664,13 @@ add_newdoc('numpy.core.multiarray', 'array', nested sequence, or if a copy is needed to satisfy any of the other requirements (`dtype`, `order`, etc.). order : {'C', 'F', 'A'}, optional - Specify the order of the array. If order is 'C' (default), then the - array will be in C-contiguous order (last-index varies the - fastest). If order is 'F', then the returned array - will be in Fortran-contiguous order (first-index varies the - fastest). If order is 'A', then the returned array may - be in any order (either C-, Fortran-contiguous, or even - discontiguous). + Specify the order of the array. If order is 'C', then the array + will be in C-contiguous order (last-index varies the fastest). + If order is 'F', then the returned array will be in + Fortran-contiguous order (first-index varies the fastest). + If order is 'A' (default), then the returned array may be + in any order (either C-, Fortran-contiguous, or even discontiguous), + unless a copy is required, in which case it will be C-contiguous. subok : bool, optional If True, then sub-classes will be passed-through, otherwise the returned array will be forced to be a base-class array (default). @@ -885,7 +885,7 @@ add_newdoc('numpy.core.multiarray', 'zeros', >>> np.zeros(5) array([ 0., 0., 0., 0., 0.]) - >>> np.zeros((5,), dtype=numpy.int) + >>> np.zeros((5,), dtype=np.int) array([0, 0, 0, 0, 0]) >>> np.zeros((2, 1)) diff --git a/numpy/core/code_generators/generate_umath.py b/numpy/core/code_generators/generate_umath.py index 9dbeb76cd..e3c4fbaa2 100644 --- a/numpy/core/code_generators/generate_umath.py +++ b/numpy/core/code_generators/generate_umath.py @@ -33,9 +33,9 @@ class TypeDescription(object): ---------- type : str Character representing the nominal type. - func_data : str or None or FullTypeDescr, optional - The string representing the expression to insert into the data array, if - any. + func_data : str or None or FullTypeDescr or FuncNameSuffix, optional + The string representing the expression to insert into the data + array, if any. in_ : str or None, optional The typecode(s) of the inputs. out : str or None, optional @@ -127,10 +127,11 @@ class Ufunc(object): import string if sys.version_info[0] < 3: - UPPER_TABLE = string.maketrans(string.ascii_lowercase, string.ascii_uppercase) + UPPER_TABLE = string.maketrans(string.ascii_lowercase, + string.ascii_uppercase) else: UPPER_TABLE = bytes.maketrans(bytes(string.ascii_lowercase, "ascii"), - bytes(string.ascii_uppercase, "ascii")) + bytes(string.ascii_uppercase, "ascii")) def english_upper(s): """ Apply English case rules to convert ASCII strings to all upper case. @@ -151,7 +152,8 @@ def english_upper(s): Examples -------- >>> from numpy.lib.utils import english_upper - >>> english_upper('ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789_') + >>> s = 'ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789_' + >>> english_upper(s) 'ABCDEFGHIJKLMNOPQRSTUVWXYZABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789_' >>> english_upper('') '' @@ -870,9 +872,9 @@ chartotype2 = {'e': 'ee_e', # 3) add function. def make_arrays(funcdict): - # functions array contains an entry for every type implemented - # NULL should be placed where PyUfunc_ style function will be filled in later - # + # functions array contains an entry for every type implemented NULL + # should be placed where PyUfunc_ style function will be filled in + # later code1list = [] code2list = [] names = sorted(funcdict.keys()) @@ -891,24 +893,25 @@ def make_arrays(funcdict): thedict = chartotype1 # one input and one output for t in uf.type_descriptions: - if t.func_data not in (None, FullTypeDescr) and not isinstance(t.func_data, FuncNameSuffix): + if (t.func_data not in (None, FullTypeDescr) and + not isinstance(t.func_data, FuncNameSuffix)): funclist.append('NULL') astype = '' if not t.astype is None: astype = '_As_%s' % thedict[t.astype] - astr = '%s_functions[%d] = PyUFunc_%s%s;' % \ - (name, k, thedict[t.type], astype) + astr = ('%s_functions[%d] = PyUFunc_%s%s;' % + (name, k, thedict[t.type], astype)) code2list.append(astr) if t.type == 'O': - astr = '%s_data[%d] = (void *) %s;' % \ - (name, k, t.func_data) + astr = ('%s_data[%d] = (void *) %s;' % + (name, k, t.func_data)) code2list.append(astr) datalist.append('(void *)NULL') elif t.type == 'P': datalist.append('(void *)"%s"' % t.func_data) else: - astr = '%s_data[%d] = (void *) %s;' % \ - (name, k, t.func_data) + astr = ('%s_data[%d] = (void *) %s;' % + (name, k, t.func_data)) code2list.append(astr) datalist.append('(void *)NULL') #datalist.append('(void *)%s' % t.func_data) @@ -916,11 +919,13 @@ def make_arrays(funcdict): elif t.func_data is FullTypeDescr: tname = english_upper(chartoname[t.type]) datalist.append('(void *)NULL') - funclist.append('%s_%s_%s_%s' % (tname, t.in_, t.out, name)) + funclist.append( + '%s_%s_%s_%s' % (tname, t.in_, t.out, name)) elif isinstance(t.func_data, FuncNameSuffix): datalist.append('(void *)NULL') tname = english_upper(chartoname[t.type]) - funclist.append('%s_%s_%s' % (tname, name, t.func_data.suffix)) + funclist.append( + '%s_%s_%s' % (tname, name, t.func_data.suffix)) else: datalist.append('(void *)NULL') tname = english_upper(chartoname[t.type]) @@ -934,11 +939,11 @@ def make_arrays(funcdict): funcnames = ', '.join(funclist) signames = ', '.join(siglist) datanames = ', '.join(datalist) - code1list.append("static PyUFuncGenericFunction %s_functions[] = { %s };" \ + code1list.append("static PyUFuncGenericFunction %s_functions[] = {%s};" % (name, funcnames)) - code1list.append("static void * %s_data[] = { %s };" \ + code1list.append("static void * %s_data[] = {%s};" % (name, datanames)) - code1list.append("static char %s_signatures[] = { %s };" \ + code1list.append("static char %s_signatures[] = {%s};" % (name, signames)) return "\n".join(code1list), "\n".join(code2list) @@ -972,8 +977,7 @@ r"""f = PyUFunc_FromFuncAndData(%s_functions, %s_data, %s_signatures, %d, name, docstring)) if uf.typereso != None: mlist.append( - r"((PyUFuncObject *)f)->type_resolver = &%s;" % - uf.typereso) + r"((PyUFuncObject *)f)->type_resolver = &%s;" % uf.typereso) mlist.append(r"""PyDict_SetItemString(dictionary, "%s", f);""" % name) mlist.append(r"""Py_DECREF(f);""") code3list.append('\n'.join(mlist)) diff --git a/numpy/core/code_generators/ufunc_docstrings.py b/numpy/core/code_generators/ufunc_docstrings.py index 328e43ca6..bfa3ad221 100644 --- a/numpy/core/code_generators/ufunc_docstrings.py +++ b/numpy/core/code_generators/ufunc_docstrings.py @@ -2310,7 +2310,7 @@ add_newdoc('numpy.core.umath', 'fmax', Returns ------- y : ndarray or scalar - The minimum of `x1` and `x2`, element-wise. Returns scalar if + The maximum of `x1` and `x2`, element-wise. Returns scalar if both `x1` and `x2` are scalars. See Also diff --git a/numpy/core/include/numpy/fenv/fenv.c b/numpy/core/include/numpy/fenv/fenv.c deleted file mode 100644 index 9a8d1be10..000000000 --- a/numpy/core/include/numpy/fenv/fenv.c +++ /dev/null @@ -1,38 +0,0 @@ -/*- - * Copyright (c) 2004 David Schultz <das@FreeBSD.ORG> - * All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * - * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - * - * $FreeBSD$ - */ - -#include <sys/types.h> -#include "fenv.h" - -const fenv_t npy__fe_dfl_env = { - 0xffff0000, - 0xffff0000, - 0xffffffff, - { 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, - 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0xff, 0xff } -}; diff --git a/numpy/core/include/numpy/fenv/fenv.h b/numpy/core/include/numpy/fenv/fenv.h deleted file mode 100644 index 79a215fc3..000000000 --- a/numpy/core/include/numpy/fenv/fenv.h +++ /dev/null @@ -1,224 +0,0 @@ -/*- - * Copyright (c) 2004 David Schultz <das@FreeBSD.ORG> - * All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * - * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - * - * $FreeBSD$ - */ - -#ifndef _FENV_H_ -#define _FENV_H_ - -#include <sys/cdefs.h> -#include <sys/types.h> - -typedef struct { - __uint32_t __control; - __uint32_t __status; - __uint32_t __tag; - char __other[16]; -} fenv_t; - -typedef __uint16_t fexcept_t; - -/* Exception flags */ -#define FE_INVALID 0x01 -#define FE_DENORMAL 0x02 -#define FE_DIVBYZERO 0x04 -#define FE_OVERFLOW 0x08 -#define FE_UNDERFLOW 0x10 -#define FE_INEXACT 0x20 -#define FE_ALL_EXCEPT (FE_DIVBYZERO | FE_DENORMAL | FE_INEXACT | \ - FE_INVALID | FE_OVERFLOW | FE_UNDERFLOW) - -/* Rounding modes */ -#define FE_TONEAREST 0x0000 -#define FE_DOWNWARD 0x0400 -#define FE_UPWARD 0x0800 -#define FE_TOWARDZERO 0x0c00 -#define _ROUND_MASK (FE_TONEAREST | FE_DOWNWARD | \ - FE_UPWARD | FE_TOWARDZERO) - -__BEGIN_DECLS - -/* Default floating-point environment */ -extern const fenv_t npy__fe_dfl_env; -#define FE_DFL_ENV (&npy__fe_dfl_env) - -#define __fldcw(__cw) __asm __volatile("fldcw %0" : : "m" (__cw)) -#define __fldenv(__env) __asm __volatile("fldenv %0" : : "m" (__env)) -#define __fnclex() __asm __volatile("fnclex") -#define __fnstenv(__env) __asm __volatile("fnstenv %0" : "=m" (*(__env))) -#define __fnstcw(__cw) __asm __volatile("fnstcw %0" : "=m" (*(__cw))) -#define __fnstsw(__sw) __asm __volatile("fnstsw %0" : "=am" (*(__sw))) -#define __fwait() __asm __volatile("fwait") - -static __inline int -feclearexcept(int __excepts) -{ - fenv_t __env; - - if (__excepts == FE_ALL_EXCEPT) { - __fnclex(); - } else { - __fnstenv(&__env); - __env.__status &= ~__excepts; - __fldenv(__env); - } - return (0); -} - -static __inline int -fegetexceptflag(fexcept_t *__flagp, int __excepts) -{ - __uint16_t __status; - - __fnstsw(&__status); - *__flagp = __status & __excepts; - return (0); -} - -static __inline int -fesetexceptflag(const fexcept_t *__flagp, int __excepts) -{ - fenv_t __env; - - __fnstenv(&__env); - __env.__status &= ~__excepts; - __env.__status |= *__flagp & __excepts; - __fldenv(__env); - return (0); -} - -static __inline int -feraiseexcept(int __excepts) -{ - fexcept_t __ex = __excepts; - - fesetexceptflag(&__ex, __excepts); - __fwait(); - return (0); -} - -static __inline int -fetestexcept(int __excepts) -{ - __uint16_t __status; - - __fnstsw(&__status); - return (__status & __excepts); -} - -static __inline int -fegetround(void) -{ - int __control; - - __fnstcw(&__control); - return (__control & _ROUND_MASK); -} - -static __inline int -fesetround(int __round) -{ - int __control; - - if (__round & ~_ROUND_MASK) - return (-1); - __fnstcw(&__control); - __control &= ~_ROUND_MASK; - __control |= __round; - __fldcw(__control); - return (0); -} - -static __inline int -fegetenv(fenv_t *__envp) -{ - int __control; - - /* - * fnstenv masks all exceptions, so we need to save and - * restore the control word to avoid this side effect. - */ - __fnstcw(&__control); - __fnstenv(__envp); - __fldcw(__control); - return (0); -} - -static __inline int -feholdexcept(fenv_t *__envp) -{ - - __fnstenv(__envp); - __fnclex(); - return (0); -} - -static __inline int -fesetenv(const fenv_t *__envp) -{ - - __fldenv(*__envp); - return (0); -} - -static __inline int -feupdateenv(const fenv_t *__envp) -{ - __uint16_t __status; - - __fnstsw(&__status); - __fldenv(*__envp); - feraiseexcept(__status & FE_ALL_EXCEPT); - return (0); -} - -#if __BSD_VISIBLE - -static __inline int -fesetmask(int __mask) -{ - int __control; - - __fnstcw(&__control); - __mask = (__control | FE_ALL_EXCEPT) & ~__mask; - __fldcw(__mask); - return (~__control & FE_ALL_EXCEPT); -} - -static __inline int -fegetmask(void) -{ - int __control; - - __fnstcw(&__control); - return (~__control & FE_ALL_EXCEPT); -} - -#endif /* __BSD_VISIBLE */ - -__END_DECLS - -#endif /* !_FENV_H_ */ diff --git a/numpy/core/include/numpy/npy_math.h b/numpy/core/include/numpy/npy_math.h index 0137e0556..8f5d8e6c7 100644 --- a/numpy/core/include/numpy/npy_math.h +++ b/numpy/core/include/numpy/npy_math.h @@ -259,7 +259,7 @@ float npy_nextafterf(float x, float y); float npy_spacingf(float x); /* - * float C99 math functions + * long double C99 math functions */ npy_longdouble npy_sinl(npy_longdouble x); @@ -425,6 +425,19 @@ npy_cdouble npy_csqrt(npy_cdouble z); npy_cdouble npy_ccos(npy_cdouble z); npy_cdouble npy_csin(npy_cdouble z); +npy_cdouble npy_ctan(npy_cdouble z); + +npy_cdouble npy_ccosh(npy_cdouble z); +npy_cdouble npy_csinh(npy_cdouble z); +npy_cdouble npy_ctanh(npy_cdouble z); + +npy_cdouble npy_cacos(npy_cdouble z); +npy_cdouble npy_casin(npy_cdouble z); +npy_cdouble npy_catan(npy_cdouble z); + +npy_cdouble npy_cacosh(npy_cdouble z); +npy_cdouble npy_casinh(npy_cdouble z); +npy_cdouble npy_catanh(npy_cdouble z); /* * Single precision complex functions @@ -440,6 +453,20 @@ npy_cfloat npy_csqrtf(npy_cfloat z); npy_cfloat npy_ccosf(npy_cfloat z); npy_cfloat npy_csinf(npy_cfloat z); +npy_cfloat npy_ctanf(npy_cfloat z); + +npy_cfloat npy_ccoshf(npy_cfloat z); +npy_cfloat npy_csinhf(npy_cfloat z); +npy_cfloat npy_ctanhf(npy_cfloat z); + +npy_cfloat npy_cacosf(npy_cfloat z); +npy_cfloat npy_casinf(npy_cfloat z); +npy_cfloat npy_catanf(npy_cfloat z); + +npy_cfloat npy_cacoshf(npy_cfloat z); +npy_cfloat npy_casinhf(npy_cfloat z); +npy_cfloat npy_catanhf(npy_cfloat z); + /* * Extended precision complex functions @@ -455,6 +482,20 @@ npy_clongdouble npy_csqrtl(npy_clongdouble z); npy_clongdouble npy_ccosl(npy_clongdouble z); npy_clongdouble npy_csinl(npy_clongdouble z); +npy_clongdouble npy_ctanl(npy_clongdouble z); + +npy_clongdouble npy_ccoshl(npy_clongdouble z); +npy_clongdouble npy_csinhl(npy_clongdouble z); +npy_clongdouble npy_ctanhl(npy_clongdouble z); + +npy_clongdouble npy_cacosl(npy_clongdouble z); +npy_clongdouble npy_casinl(npy_clongdouble z); +npy_clongdouble npy_catanl(npy_clongdouble z); + +npy_clongdouble npy_cacoshl(npy_clongdouble z); +npy_clongdouble npy_casinhl(npy_clongdouble z); +npy_clongdouble npy_catanhl(npy_clongdouble z); + /* * Functions that set the floating point error diff --git a/numpy/core/numeric.py b/numpy/core/numeric.py index a464562c4..7a0fa4b62 100644 --- a/numpy/core/numeric.py +++ b/numpy/core/numeric.py @@ -1024,9 +1024,8 @@ def outer(a, b, out=None): Second input vector. Input is flattened if not already 1-dimensional. out : (M, N) ndarray, optional - A location where the result is stored - .. versionadded:: 1.9.0 + A location where the result is stored Returns ------- @@ -2206,42 +2205,7 @@ def identity(n, dtype=None): from numpy import eye return eye(n, dtype=dtype) -def _allclose_points(a, b, rtol=1.e-5, atol=1.e-8): - """ - This is the point-wise inner calculation of 'allclose', which is subtly - different from 'isclose'. Provided as a comparison routine for use in - testing.assert_allclose. - See those routines for further details. - - """ - x = array(a, copy=False, ndmin=1) - y = array(b, copy=False, ndmin=1) - - # make sure y is an inexact type to avoid abs(MIN_INT); will cause - # casting of x later. - dtype = multiarray.result_type(y, 1.) - y = array(y, dtype=dtype, copy=False) - - xinf = isinf(x) - yinf = isinf(y) - if any(xinf) or any(yinf): - # Check that x and y have inf's only in the same positions - if not all(xinf == yinf): - return False - # Check that sign of inf's in x and y is the same - if not all(x[xinf] == y[xinf]): - return False - - x = x[~xinf] - y = y[~xinf] - - # ignore invalid fpe's - with errstate(invalid='ignore'): - r = less_equal(abs(x - y), atol + rtol * abs(y)) - - return r - -def allclose(a, b, rtol=1.e-5, atol=1.e-8): +def allclose(a, b, rtol=1.e-5, atol=1.e-8, equal_nan=False): """ Returns True if two arrays are element-wise equal within a tolerance. @@ -2262,6 +2226,10 @@ def allclose(a, b, rtol=1.e-5, atol=1.e-8): The relative tolerance parameter (see Notes). atol : float The absolute tolerance parameter (see Notes). + equal_nan : bool + .. versionadded:: 1.10.0 + Whether to compare NaN's as equal. If True, NaN's in `a` will be + considered equal to NaN's in `b` in the output array. Returns ------- @@ -2294,9 +2262,11 @@ def allclose(a, b, rtol=1.e-5, atol=1.e-8): False >>> np.allclose([1.0, np.nan], [1.0, np.nan]) False + >>> np.allclose([1.0, np.nan], [1.0, np.nan], equal_nan=True) + True """ - return all(_allclose_points(a, b, rtol=rtol, atol=atol)) + return all(isclose(a, b, rtol=rtol, atol=atol, equal_nan=equal_nan)) def isclose(a, b, rtol=1.e-5, atol=1.e-8, equal_nan=False): """ @@ -2366,6 +2336,13 @@ def isclose(a, b, rtol=1.e-5, atol=1.e-8, equal_nan=False): x = array(a, copy=False, subok=True, ndmin=1) y = array(b, copy=False, subok=True, ndmin=1) + + # Make sure y is an inexact type to avoid bad behavior on abs(MIN_INT). + # This will cause casting of x later. Also, make sure to allow subclasses + # (e.g., for numpy.ma). + dt = multiarray.result_type(y, 1.) + y = array(y, dtype=dt, copy=False, subok=True) + xfin = isfinite(x) yfin = isfinite(y) if all(xfin) and all(yfin): diff --git a/numpy/core/records.py b/numpy/core/records.py index 23680711f..243645436 100644 --- a/numpy/core/records.py +++ b/numpy/core/records.py @@ -40,7 +40,6 @@ import sys import os from . import numeric as sb -from .defchararray import chararray from . import numerictypes as nt from numpy.compat import isfileobj, bytes, long @@ -215,6 +214,12 @@ class format_parser: class record(nt.void): """A data-type scalar that allows field access as attribute lookup. """ + + # manually set name and module so that this class's type shows up + # as numpy.record when printed + __name__ = 'record' + __module__ = 'numpy' + def __repr__(self): return self.__str__() @@ -232,17 +237,15 @@ class record(nt.void): res = fielddict.get(attr, None) if res: obj = self.getfield(*res[:2]) - # if it has fields return a recarray, - # if it's a string ('SU') return a chararray + # if it has fields return a record, # otherwise return the object try: dt = obj.dtype except AttributeError: + #happens if field is Object type return obj if dt.fields: - return obj.view(obj.__class__) - if dt.char in 'SU': - return obj.view(chararray) + return obj.view((record, obj.dtype.descr)) return obj else: raise AttributeError("'record' object has no " @@ -388,6 +391,12 @@ class recarray(ndarray): dtype=[('x', '<i4'), ('y', '<f8'), ('z', '<i4')]) """ + + # manually set name and module so that this class's type shows + # up as "numpy.recarray" when printed + __name__ = 'recarray' + __module__ = 'numpy' + def __new__(subtype, shape, dtype=None, buf=None, offset=0, strides=None, formats=None, names=None, titles=None, byteorder=None, aligned=False, order='C'): @@ -406,29 +415,37 @@ class recarray(ndarray): return self def __getattribute__(self, attr): + # See if ndarray has this attr, and return it if so. (note that this + # means a field with the same name as an ndarray attr cannot be + # accessed by attribute). try: return object.__getattribute__(self, attr) except AttributeError: # attr must be a fieldname pass + + # look for a field with this name fielddict = ndarray.__getattribute__(self, 'dtype').fields try: res = fielddict[attr][:2] except (TypeError, KeyError): - raise AttributeError("record array has no attribute %s" % attr) + raise AttributeError("recarray has no attribute %s" % attr) obj = self.getfield(*res) - # if it has fields return a recarray, otherwise return - # normal array - if obj.dtype.fields: - return obj - if obj.dtype.char in 'SU': - return obj.view(chararray) - return obj.view(ndarray) -# Save the dictionary -# If the attr is a field name and not in the saved dictionary -# Undo any "setting" of the attribute and do a setfield -# Thus, you can't create attributes on-the-fly that are field names. + # At this point obj will always be a recarray, since (see + # PyArray_GetField) the type of obj is inherited. Next, if obj.dtype is + # non-structured, convert it to an ndarray. If obj is structured leave + # it as a recarray, but make sure to convert to the same dtype.type (eg + # to preserve numpy.record type if present), since nested structured + # fields do not inherit type. + if obj.dtype.fields: + return obj.view(dtype=(self.dtype.type, obj.dtype.descr)) + else: + return obj.view(ndarray) + # Save the dictionary. + # If the attr is a field name and not in the saved dictionary + # Undo any "setting" of the attribute and do a setfield + # Thus, you can't create attributes on-the-fly that are field names. def __setattr__(self, attr, val): newattr = attr not in self.__dict__ try: @@ -456,13 +473,40 @@ class recarray(ndarray): def __getitem__(self, indx): obj = ndarray.__getitem__(self, indx) - if (isinstance(obj, ndarray) and obj.dtype.isbuiltin): - return obj.view(ndarray) - return obj - def __repr__(self) : - ret = ndarray.__repr__(self) - return ret.replace("recarray", "rec.array", 1) + # copy behavior of getattr, except that here + # we might also be returning a single element + if isinstance(obj, ndarray): + if obj.dtype.fields: + return obj.view(dtype=(self.dtype.type, obj.dtype.descr)) + else: + return obj.view(type=ndarray) + else: + # return a single element + return obj + + def __repr__(self): + # get data/shape string. logic taken from numeric.array_repr + if self.size > 0 or self.shape==(0,): + lst = sb.array2string(self, separator=', ') + else: + # show zero-length shape unless it is (0,) + lst = "[], shape=%s" % (repr(self.shape),) + + if self.dtype.type is record: + # If this is a full record array (has numpy.record dtype), + # represent it using the rec.array function. Since rec.array + # converts dtype to a numpy.record for us, use only dtype.descr, + # not repr(dtype). + lf = '\n'+' '*len("rec.array(") + return ('rec.array(%s, %sdtype=%s)' % + (lst, lf, repr(self.dtype.descr))) + else: + # otherwise represent it using np.array plus a view + # (There is currently (v1.10) no other easy way to create it) + lf = '\n'+' '*len("array(") + return ('array(%s, %sdtype=%s).view(numpy.recarray)' % + (lst, lf, str(self.dtype))) def field(self, attr, val=None): if isinstance(attr, int): @@ -477,8 +521,6 @@ class recarray(ndarray): obj = self.getfield(*res) if obj.dtype.fields: return obj - if obj.dtype.char in 'SU': - return obj.view(chararray) return obj.view(ndarray) else: return self.setfield(val, *res) @@ -589,7 +631,7 @@ def fromrecords(recList, dtype=None, shape=None, formats=None, names=None, >>> r.col1 array([456, 2]) >>> r.col2 - chararray(['dbe', 'de'], + array(['dbe', 'de'], dtype='|S3') >>> import pickle >>> print pickle.loads(pickle.dumps(r)) diff --git a/numpy/core/setup.py b/numpy/core/setup.py index 75d64d81b..7f0649158 100644 --- a/numpy/core/setup.py +++ b/numpy/core/setup.py @@ -630,15 +630,9 @@ def configuration(parent_package='',top_path=None): deps = [join('src', 'npymath', '_signbit.c'), join('include', 'numpy', '*object.h'), - 'include/numpy/fenv/fenv.c', - 'include/numpy/fenv/fenv.h', join(codegen_dir, 'genapi.py'), ] - # Don't install fenv unless we need them. - if sys.platform == 'cygwin': - config.add_data_dir('include/numpy/fenv') - ####################################################################### # dummy module # ####################################################################### diff --git a/numpy/core/setup_common.py b/numpy/core/setup_common.py index 0b18bc6c6..6abbe5ff2 100644 --- a/numpy/core/setup_common.py +++ b/numpy/core/setup_common.py @@ -107,6 +107,7 @@ OPTIONAL_HEADERS = [ # sse headers only enabled automatically on amd64/x32 builds "xmmintrin.h", # SSE "emmintrin.h", # SSE2 + "features.h", # for glibc version linux ] # optional gcc compiler builtins and their call arguments and optional a @@ -138,23 +139,29 @@ OPTIONAL_FUNCTION_ATTRIBUTES = [('__attribute__((optimize("unroll-loops")))', OPTIONAL_VARIABLE_ATTRIBUTES = ["__thread", "__declspec(thread)"] # Subset of OPTIONAL_STDFUNCS which may alreay have HAVE_* defined by Python.h -OPTIONAL_STDFUNCS_MAYBE = ["expm1", "log1p", "acosh", "atanh", "asinh", "hypot", - "copysign", "ftello", "fseeko"] +OPTIONAL_STDFUNCS_MAYBE = [ + "expm1", "log1p", "acosh", "atanh", "asinh", "hypot", "copysign", + "ftello", "fseeko" + ] # C99 functions: float and long double versions -C99_FUNCS = ["sin", "cos", "tan", "sinh", "cosh", "tanh", "fabs", "floor", - "ceil", "rint", "trunc", "sqrt", "log10", "log", "log1p", "exp", - "expm1", "asin", "acos", "atan", "asinh", "acosh", "atanh", - "hypot", "atan2", "pow", "fmod", "modf", 'frexp', 'ldexp', - "exp2", "log2", "copysign", "nextafter", "cbrt"] - +C99_FUNCS = [ + "sin", "cos", "tan", "sinh", "cosh", "tanh", "fabs", "floor", "ceil", + "rint", "trunc", "sqrt", "log10", "log", "log1p", "exp", "expm1", + "asin", "acos", "atan", "asinh", "acosh", "atanh", "hypot", "atan2", + "pow", "fmod", "modf", 'frexp', 'ldexp', "exp2", "log2", "copysign", + "nextafter", "cbrt" + ] C99_FUNCS_SINGLE = [f + 'f' for f in C99_FUNCS] C99_FUNCS_EXTENDED = [f + 'l' for f in C99_FUNCS] - -C99_COMPLEX_TYPES = ['complex double', 'complex float', 'complex long double'] - -C99_COMPLEX_FUNCS = ['creal', 'cimag', 'cabs', 'carg', 'cexp', 'csqrt', 'clog', - 'ccos', 'csin', 'cpow'] +C99_COMPLEX_TYPES = [ + 'complex double', 'complex float', 'complex long double' + ] +C99_COMPLEX_FUNCS = [ + "cabs", "cacos", "cacosh", "carg", "casin", "casinh", "catan", + "catanh", "ccos", "ccosh", "cexp", "cimag", "clog", "conj", "cpow", + "cproj", "creal", "csin", "csinh", "csqrt", "ctan", "ctanh" + ] def fname2def(name): return "HAVE_%s" % name.upper() diff --git a/numpy/core/src/multiarray/cblasfuncs.c b/numpy/core/src/multiarray/cblasfuncs.c index f8f508ed0..9c7c26725 100644 --- a/numpy/core/src/multiarray/cblasfuncs.c +++ b/numpy/core/src/multiarray/cblasfuncs.c @@ -685,7 +685,7 @@ cblas_innerproduct(int typenum, PyArrayObject *ap1, PyArrayObject *ap2) int j, l, lda, ldb; int nd; double prior1, prior2; - PyArrayObject *ret; + PyArrayObject *ret = NULL; npy_intp dimensions[NPY_MAXDIMS]; PyTypeObject *subtype; diff --git a/numpy/core/src/multiarray/descriptor.c b/numpy/core/src/multiarray/descriptor.c index 9cf66020d..0993190b7 100644 --- a/numpy/core/src/multiarray/descriptor.c +++ b/numpy/core/src/multiarray/descriptor.c @@ -3159,6 +3159,8 @@ arraydescr_struct_dict_str(PyArray_Descr *dtype, int includealignedflag) static PyObject * arraydescr_struct_str(PyArray_Descr *dtype, int includealignflag) { + PyObject *sub; + /* * The list str representation can't include the 'align=' flag, * so if it is requested and the struct has the aligned flag set, @@ -3166,10 +3168,52 @@ arraydescr_struct_str(PyArray_Descr *dtype, int includealignflag) */ if (!(includealignflag && (dtype->flags&NPY_ALIGNED_STRUCT)) && is_dtype_struct_simple_unaligned_layout(dtype)) { - return arraydescr_struct_list_str(dtype); + sub = arraydescr_struct_list_str(dtype); + } + else { + sub = arraydescr_struct_dict_str(dtype, includealignflag); + } + + /* If the data type has a non-void (subclassed) type, show it */ + if (dtype->type_num == NPY_VOID && dtype->typeobj != &PyVoidArrType_Type) { + /* + * Note: We cannot get the type name from dtype->typeobj->tp_name + * because its value depends on whether the type is dynamically or + * statically allocated. Instead use __name__ and __module__. + * See https://docs.python.org/2/c-api/typeobj.html. + */ + + PyObject *str_name, *namestr, *str_module, *modulestr, *ret; + + str_name = PyUString_FromString("__name__"); + namestr = PyObject_GetAttr((PyObject*)(dtype->typeobj), str_name); + Py_DECREF(str_name); + + if (namestr == NULL) { + /* this should never happen since types always have __name__ */ + PyErr_Format(PyExc_RuntimeError, + "dtype does not have a __name__ attribute"); + return NULL; + } + + str_module = PyUString_FromString("__module__"); + modulestr = PyObject_GetAttr((PyObject*)(dtype->typeobj), str_module); + Py_DECREF(str_module); + + ret = PyUString_FromString("("); + if (modulestr != NULL) { + /* Note: if modulestr == NULL, the type is unpicklable */ + PyUString_ConcatAndDel(&ret, modulestr); + PyUString_ConcatAndDel(&ret, PyUString_FromString(".")); + } + PyUString_ConcatAndDel(&ret, namestr); + PyUString_ConcatAndDel(&ret, PyUString_FromString(", ")); + PyUString_ConcatAndDel(&ret, sub); + PyUString_ConcatAndDel(&ret, PyUString_FromString(")")); + return ret; } else { - return arraydescr_struct_dict_str(dtype, includealignflag); + return sub; } } diff --git a/numpy/core/src/multiarray/scalartypes.c.src b/numpy/core/src/multiarray/scalartypes.c.src index 4099326f1..252e5b726 100644 --- a/numpy/core/src/multiarray/scalartypes.c.src +++ b/numpy/core/src/multiarray/scalartypes.c.src @@ -114,12 +114,17 @@ gentype_alloc(PyTypeObject *type, Py_ssize_t nitems) const size_t size = _PyObject_VAR_SIZE(type, nitems + 1); obj = (PyObject *)PyObject_Malloc(size); + /* + * Fixme. Need to check for no memory. + * If we don't need to zero memory, we could use + * PyObject_{New, NewVar} for this whole function. + */ memset(obj, 0, size); if (type->tp_itemsize == 0) { - PyObject_INIT(obj, type); + PyObject_Init(obj, type); } else { - (void) PyObject_INIT_VAR((PyVarObject *)obj, type, nitems); + (void) PyObject_InitVar((PyVarObject *)obj, type, nitems); } return obj; } diff --git a/numpy/core/src/multiarray/shape.c b/numpy/core/src/multiarray/shape.c index f1e81ff6b..03bfc6a7a 100644 --- a/numpy/core/src/multiarray/shape.c +++ b/numpy/core/src/multiarray/shape.c @@ -948,7 +948,7 @@ PyArray_Ravel(PyArrayObject *arr, NPY_ORDER order) /* For KEEPORDER, check if we can make a flattened view */ else { npy_stride_sort_item strideperm[NPY_MAXDIMS]; - npy_intp stride, base_stride = NPY_MIN_INTP; + npy_intp stride = 0, base_stride = NPY_MIN_INTP; int i, ndim = PyArray_NDIM(arr); PyArray_CreateSortedStridePerm(PyArray_NDIM(arr), diff --git a/numpy/core/src/npymath/ieee754.c.src b/numpy/core/src/npymath/ieee754.c.src index b8f6889e1..0370ea6c7 100644 --- a/numpy/core/src/npymath/ieee754.c.src +++ b/numpy/core/src/npymath/ieee754.c.src @@ -615,13 +615,7 @@ void npy_set_floatstatus_invalid(void) #elif defined(__GLIBC__) || defined(__APPLE__) || \ defined(__CYGWIN__) || defined(__MINGW32__) || \ (defined(__FreeBSD__) && (__FreeBSD_version >= 502114)) - -# if defined(__GLIBC__) || defined(__APPLE__) || \ - defined(__MINGW32__) || defined(__FreeBSD__) # include <fenv.h> -# elif defined(__CYGWIN__) -# include "numpy/fenv/fenv.h" -# endif int npy_get_floatstatus(void) { diff --git a/numpy/core/src/npymath/npy_math_complex.c.src b/numpy/core/src/npymath/npy_math_complex.c.src index 9cbfd32ae..94cb7e29d 100644 --- a/numpy/core/src/npymath/npy_math_complex.c.src +++ b/numpy/core/src/npymath/npy_math_complex.c.src @@ -1,10 +1,13 @@ /* + * vim: syntax=c + * * Implement some C99-compatible complex math functions * - * Most of the code is taken from the msun library in FreeBSD (HEAD @ 30th June - * 2009), under the following license: + * Most of the code is taken from the msun library in FreeBSD (HEAD @ 4th + * October 2013), under the following license: * - * Copyright (c) 2007 David Schultz <das@FreeBSD.ORG> + * Copyright (c) 2007, 2011 David Schultz <das@FreeBSD.ORG> + * Copyright (c) 2012 Stephen Montgomery-Smith <stephen@FreeBSD.ORG> * All rights reserved. * * Redistribution and use in source and binary forms, with or without @@ -31,9 +34,12 @@ #include "npy_math_common.h" #include "npy_math_private.h" -/*========================================================== - * Custom implementation of missing complex C99 functions - *=========================================================*/ + +#define raise_inexact() do { volatile npy_float junk = 1 + tiny; } while(0) + + +static npy_float tiny = 3.9443045e-31f; + /**begin repeat * #type = npy_float, npy_double, npy_longdouble# @@ -41,23 +47,171 @@ * #c = f, , l# * #C = F, , L# * #TMAX = FLT_MAX, DBL_MAX, LDBL_MAX# + * #TMIN = FLT_MIN, DBL_MIN, LDBL_MIN# + * #TMANT_DIG = FLT_MANT_DIG, DBL_MANT_DIG, LDBL_MANT_DIG# + * #TEPS = FLT_EPSILON, DBL_EPSILON, LDBL_EPSILON# + * #precision = 1, 2, 3# */ + +/*========================================================== + * Constants + *=========================================================*/ +static const @ctype@ c_1@c@ = {1.0@C@, 0.0}; +static const @ctype@ c_half@c@ = {0.5@C@, 0.0}; +static const @ctype@ c_i@c@ = {0.0, 1.0@C@}; +static const @ctype@ c_ihalf@c@ = {0.0, 0.5@C@}; + +/*========================================================== + * Helper functions + * + * These are necessary because we do not count on using a + * C99 compiler. + *=========================================================*/ +static NPY_INLINE +@ctype@ +cadd@c@(@ctype@ a, @ctype@ b) +{ + return npy_cpack@c@(npy_creal@c@(a) + npy_creal@c@(b), + npy_cimag@c@(a) + npy_cimag@c@(b)); +} + +static NPY_INLINE +@ctype@ +csub@c@(@ctype@ a, @ctype@ b) +{ + return npy_cpack@c@(npy_creal@c@(a) - npy_creal@c@(b), + npy_cimag@c@(a) - npy_cimag@c@(b)); +} + +static NPY_INLINE +@ctype@ +cmul@c@(@ctype@ a, @ctype@ b) +{ + @type@ ar, ai, br, bi; + ar = npy_creal@c@(a); + ai = npy_cimag@c@(a); + br = npy_creal@c@(b); + bi = npy_cimag@c@(b); + return npy_cpack@c@(ar*br - ai*bi, ar*bi + ai*br); +} + +static NPY_INLINE +@ctype@ +cdiv@c@(@ctype@ a, @ctype@ b) +{ + @type@ ar, ai, br, bi, abs_br, abs_bi; + ar = npy_creal@c@(a); + ai = npy_cimag@c@(a); + br = npy_creal@c@(b); + bi = npy_cimag@c@(b); + abs_br = npy_fabs@c@(br); + abs_bi = npy_fabs@c@(bi); + + if (abs_br >= abs_bi) { + if (abs_br == 0 && abs_bi == 0) { + /* divide by zeros should yield a complex inf or nan */ + return npy_cpack@c@(ar/abs_br, ai/abs_bi); + } + else { + @type@ rat = bi/br; + @type@ scl = 1.0@C@/(br+bi*rat); + return npy_cpack@c@((ar + ai*rat)*scl, (ai - ar*rat)*scl); + } + } + else { + @type@ rat = br/bi; + @type@ scl = 1.0@C@/(bi + br*rat); + return npy_cpack@c@((ar*rat + ai)*scl, (ai*rat - ar)*scl); + } +} + +static NPY_INLINE +@ctype@ +cneg@c@(@ctype@ a) +{ + return npy_cpack@c@(-npy_creal@c@(a), -npy_cimag@c@(a)); +} + +static NPY_INLINE +@ctype@ +cmuli@c@(@ctype@ a) +{ + return npy_cpack@c@(-npy_cimag@c@(a), npy_creal@c@(a)); +} + +/*========================================================== + * Custom implementation of missing complex C99 functions + *=========================================================*/ + #ifndef HAVE_CABS@C@ -@type@ npy_cabs@c@(@ctype@ z) +@type@ +npy_cabs@c@(@ctype@ z) { return npy_hypot@c@(npy_creal@c@(z), npy_cimag@c@(z)); } #endif #ifndef HAVE_CARG@C@ -@type@ npy_carg@c@(@ctype@ z) +@type@ +npy_carg@c@(@ctype@ z) { return npy_atan2@c@(npy_cimag@c@(z), npy_creal@c@(z)); } #endif +/* + * cexp and (ccos, csin)h functions need to calculate exp scaled by another + * number. This can be difficult if exp(x) overflows. By doing this way, we + * don't risk overflowing exp. This likely raises floating-point exceptions, + * if we decide that we care. + * + * This is only useful over a limited range, (see below) an expects that the + * input values are in this range. + * + * This is based on the technique used in FreeBSD's __frexp_exp and + * __ldexp_(c)exp functions by David Schultz. + * + * SCALED_CEXP_LOWER = log(FLT_MAX) + * SCALED_CEXP_UPPER = log(2) + log(FLT_MAX) - log(FLT_TRUE_MIN), + * where FLT_TRUE_MIN is the smallest possible subnormal number. + */ + +#define SCALED_CEXP_LOWERF 88.722839f +#define SCALED_CEXP_UPPERF 192.69492f +#define SCALED_CEXP_LOWER 710.47586007394386 +#define SCALED_CEXP_UPPER 1454.9159319953251 +#define SCALED_CEXP_LOWERL 11357.216553474703895L +#define SCALED_CEXP_UPPERL 22756.021937783004509L + +static +@ctype@ +_npy_scaled_cexp@c@(@type@ x, @type@ y, npy_int expt) +{ +#if @precision@ == 1 + const npy_int k = 235; +#endif +#if @precision@ == 2 + const npy_int k = 1799; +#endif +#if @precision@ == 3 + const npy_int k = 19547; +#endif + const @type@ kln2 = k * NPY_LOGE2@c@; + @type@ mant, mantcos, mantsin; + npy_int ex, excos, exsin; + + mant = npy_frexp@c@(npy_exp@c@(x - kln2), &ex); + mantcos = npy_frexp@c@(npy_cos@c@(y), &excos); + mantsin = npy_frexp@c@(npy_sin@c@(y), &exsin); + + expt += ex + k; + return npy_cpack@c@( npy_ldexp@c@(mant * mantcos, expt + excos), + npy_ldexp@c@(mant * mantsin, expt + exsin)); +} + #ifndef HAVE_CEXP@C@ -@ctype@ npy_cexp@c@(@ctype@ z) +@ctype@ +npy_cexp@c@(@ctype@ z) { @type@ x, c, s; @type@ r, i; @@ -67,46 +221,60 @@ i = npy_cimag@c@(z); if (npy_isfinite(r)) { - x = npy_exp@c@(r); + if (r >= SCALED_CEXP_LOWER@C@ && r <= SCALED_CEXP_UPPER@C@) { + ret = _npy_scaled_cexp@c@(r, i, 0); + } + else { + x = npy_exp@c@(r); - c = npy_cos@c@(i); - s = npy_sin@c@(i); + c = npy_cos@c@(i); + s = npy_sin@c@(i); - if (npy_isfinite(i)) { - ret = npy_cpack@c@(x * c, x * s); - } else { - ret = npy_cpack@c@(NPY_NAN, npy_copysign@c@(NPY_NAN, i)); + if (npy_isfinite(i)) { + ret = npy_cpack@c@(x * c, x * s); + } + else { + ret = npy_cpack@c@(NPY_NAN@C@, npy_copysign@c@(NPY_NAN@C@, i)); + } } - } else if (npy_isnan(r)) { + } + else if (npy_isnan(r)) { /* r is nan */ if (i == 0) { - ret = npy_cpack@c@(r, 0); - } else { - ret = npy_cpack@c@(r, npy_copysign@c@(NPY_NAN, i)); + ret = z; + } + else { + ret = npy_cpack@c@(r, npy_copysign@c@(NPY_NAN@C@, i)); } - } else { + } + else { /* r is +- inf */ if (r > 0) { if (i == 0) { ret = npy_cpack@c@(r, i); - } else if (npy_isfinite(i)) { + } + else if (npy_isfinite(i)) { c = npy_cos@c@(i); s = npy_sin@c@(i); ret = npy_cpack@c@(r * c, r * s); - } else { + } + else { /* x = +inf, y = +-inf | nan */ - ret = npy_cpack@c@(r, NPY_NAN); + npy_set_floatstatus_invalid(); + ret = npy_cpack@c@(r, NPY_NAN@C@); } - } else { + } + else { if (npy_isfinite(i)) { x = npy_exp@c@(r); c = npy_cos@c@(i); s = npy_sin@c@(i); ret = npy_cpack@c@(x * c, x * s); - } else { + } + else { /* x = -inf, y = nan | +i inf */ ret = npy_cpack@c@(0, 0); } @@ -118,9 +286,72 @@ #endif #ifndef HAVE_CLOG@C@ -@ctype@ npy_clog@c@(@ctype@ z) +/* algorithm from cpython, rev. d86f5686cef9 + * + * 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. +*/ +@ctype@ +npy_clog@c@(@ctype@ z) { - return npy_cpack@c@(npy_log@c@ (npy_cabs@c@ (z)), npy_carg@c@ (z)); + @type@ ax = npy_fabs@c@(npy_creal@c@(z)); + @type@ ay = npy_fabs@c@(npy_cimag@c@(z)); + @type@ rr, ri; + + if (ax > @TMAX@/4 || ay > @TMAX@/4) { + rr = npy_log@c@(npy_hypot@c@(ax/2, ay/2)) + NPY_LOGE2@c@; + } + else if (ax < @TMIN@ && ay < @TMIN@) { + if (ax > 0 || ay > 0) { + /* catch cases where hypot(ax, ay) is subnormal */ + rr = npy_log@c@(npy_hypot@c@(npy_ldexp@c@(ax, @TMANT_DIG@), + npy_ldexp@c@(ay, @TMANT_DIG@))) - @TMANT_DIG@*NPY_LOGE2@c@; + } + else { + /* log(+/-0 +/- 0i) */ + /* raise divide-by-zero floating point exception */ + rr = -1.0@c@ / npy_creal@c@(z); + rr = npy_copysign@c@(rr, -1); + ri = npy_carg@c@(z); + return npy_cpack@c@(rr, ri); + } + } + else { + @type@ h = npy_hypot@c@(ax, ay); + if (0.71 <= h && h <= 1.73) { + @type@ am = ax > ay ? ax : ay; /* max(ax, ay) */ + @type@ an = ax > ay ? ay : ax; /* min(ax, ay) */ + rr = npy_log1p@c@((am-1)*(am+1)+an*an)/2; + } + else { + rr = npy_log@c@(h); + } + } + ri = npy_carg@c@(z); + + return npy_cpack@c@(rr, ri); } #endif @@ -129,7 +360,8 @@ /* We risk spurious overflow for components >= DBL_MAX / (1 + sqrt(2)). */ #define THRESH (@TMAX@ / (1 + NPY_SQRT2@c@)) -@ctype@ npy_csqrt@c@(@ctype@ z) +@ctype@ +npy_csqrt@c@(@ctype@ z) { @ctype@ result; @type@ a, b; @@ -140,10 +372,12 @@ b = npy_cimag@c@(z); /* Handle special cases. */ - if (a == 0 && b == 0) + if (a == 0 && b == 0) { return (npy_cpack@c@(0, b)); - if (npy_isinf(b)) - return (npy_cpack@c@(NPY_INFINITY, b)); + } + if (npy_isinf(b)) { + return (npy_cpack@c@(NPY_INFINITY@C@, b)); + } if (npy_isnan(a)) { t = (b - b) / (b - b); /* raise invalid if b is not a NaN */ return (npy_cpack@c@(a, t)); /* return NaN + NaN i */ @@ -155,10 +389,12 @@ * csqrt(-inf + NaN i) = NaN +- inf i * csqrt(-inf + y i) = 0 + inf i */ - if (npy_signbit(a)) + if (npy_signbit(a)) { return (npy_cpack@c@(npy_fabs@c@(b - b), npy_copysign@c@(a, b))); - else + } + else { return (npy_cpack@c@(a, npy_copysign@c@(b - b, b))); + } } /* * The remaining special case (b is NaN) is handled just fine by @@ -170,61 +406,1333 @@ a *= 0.25; b *= 0.25; scale = 1; - } else { + } + else { scale = 0; } /* Algorithm 312, CACM vol 10, Oct 1967. */ if (a >= 0) { - t = npy_sqrt@c@((a + npy_hypot@c@(a, b)) * 0.5); + t = npy_sqrt@c@((a + npy_hypot@c@(a, b)) * 0.5@c@); result = npy_cpack@c@(t, b / (2 * t)); - } else { - t = npy_sqrt@c@((-a + npy_hypot@c@(a, b)) * 0.5); + } + else { + t = npy_sqrt@c@((-a + npy_hypot@c@(a, b)) * 0.5@c@); result = npy_cpack@c@(npy_fabs@c@(b) / (2 * t), npy_copysign@c@(t, b)); } /* Rescale. */ - if (scale) + if (scale) { return (npy_cpack@c@(npy_creal@c@(result) * 2, npy_cimag@c@(result))); - else + } + else { return (result); + } } #undef THRESH #endif -#ifndef HAVE_CPOW@C@ -@ctype@ npy_cpow@c@ (@ctype@ x, @ctype@ y) +/* + * Always use this function because of the multiplication for small + * integer powers, but in the body use cpow if it is available. + */ + +/* private function for use in npy_pow{f, ,l} */ +#ifdef HAVE_CPOW@C@ +static @ctype@ +sys_cpow@c@(@ctype@ x, @ctype@ y) { - @ctype@ b; - @type@ br, bi, yr, yi; + __@ctype@_to_c99_cast xcast; + __@ctype@_to_c99_cast ycast; + __@ctype@_to_c99_cast ret; + xcast.npy_z = x; + ycast.npy_z = y; + ret.c99_z = cpow@c@(xcast.c99_z, ycast.c99_z); + return ret.npy_z; +} +#endif - yr = npy_creal@c@(y); - yi = npy_cimag@c@(y); - b = npy_clog@c@(x); - br = npy_creal@c@(b); - bi = npy_cimag@c@(b); - return npy_cexp@c@(npy_cpack@c@(br * yr - bi * yi, br * yi + bi * yr)); -} +@ctype@ +npy_cpow@c@ (@ctype@ a, @ctype@ b) +{ + npy_intp n; + @type@ ar = npy_creal@c@(a); + @type@ br = npy_creal@c@(b); + @type@ ai = npy_cimag@c@(a); + @type@ bi = npy_cimag@c@(b); + @ctype@ r; + + if (br == 0. && bi == 0.) { + return npy_cpack@c@(1., 0.); + } + if (ar == 0. && ai == 0.) { + if (br > 0 && bi == 0) { + return npy_cpack@c@(0., 0.); + } + else { + volatile @type@ tmp = NPY_INFINITY@C@; + /* + * NB: there are four complex zeros; c0 = (+-0, +-0), so that + * unlike for reals, c0**p, with `p` negative is in general + * ill-defined. + * + * c0**z with z complex is also ill-defined. + */ + r = npy_cpack@c@(NPY_NAN@C@, NPY_NAN@C@); + + /* Raise invalid */ + tmp -= NPY_INFINITY@C@; + ar = tmp; + return r; + } + } + if (bi == 0 && (n=(npy_intp)br) == br) { + if (n == 1) { + /* unroll: handle inf better */ + return npy_cpack@c@(ar, ai); + } + else if (n == 2) { + /* unroll: handle inf better */ + return cmul@c@(a, a); + } + else if (n == 3) { + /* unroll: handle inf better */ + return cmul@c@(a, cmul@c@(a, a)); + } + else if (n > -100 && n < 100) { + @ctype@ p, aa; + npy_intp mask = 1; + if (n < 0) { + n = -n; + } + aa = c_1@c@; + p = npy_cpack@c@(ar, ai); + while (1) { + if (n & mask) { + aa = cmul@c@(aa,p); + } + mask <<= 1; + if (n < mask || mask <= 0) { + break; + } + p = cmul@c@(p,p); + } + r = npy_cpack@c@(npy_creal@c@(aa), npy_cimag@c@(aa)); + if (br < 0) { + r = cdiv@c@(c_1@c@, r); + } + return r; + } + } + +#ifdef HAVE_CPOW@C@ + return sys_cpow@c@(a, b); + +#else + { + @ctype@ loga = npy_clog@c@(a); + + ar = npy_creal@c@(loga); + ai = npy_cimag@c@(loga); + return npy_cexp@c@(npy_cpack@c@(ar*br - ai*bi, ar*bi + ai*br)); + } + #endif +} + #ifndef HAVE_CCOS@C@ -@ctype@ npy_ccos@c@(@ctype@ z) +@ctype@ +npy_ccos@c@(@ctype@ z) { - @type@ x, y; + /* ccos(z) = ccosh(I * z) */ + return npy_ccosh@c@(npy_cpack@c@(-npy_cimag@c@(z), npy_creal@c@(z))); +} +#endif + +#ifndef HAVE_CSIN@C@ +@ctype@ +npy_csin@c@(@ctype@ z) +{ + /* csin(z) = -I * csinh(I * z) */ + z = npy_csinh@c@(npy_cpack@c@(-npy_cimag@c@(z), npy_creal@c@(z))); + return npy_cpack@c@(npy_cimag@c@(z), -npy_creal@c@(z)); +} +#endif + +#ifndef HAVE_CTAN@C@ +@ctype@ +npy_ctan@c@(@ctype@ z) +{ + /* ctan(z) = -I * ctanh(I * z) */ + z = npy_ctanh@c@(npy_cpack@c@(-npy_cimag@c@(z), npy_creal@c@(z))); + return (npy_cpack@c@(npy_cimag@c@(z), -npy_creal@c@(z))); +} +#endif + +#ifndef HAVE_CCOSH@C@ +/* + * Taken from the msun library in FreeBSD, rev 226599. + * + * Hyperbolic cosine of a complex argument z = x + i y. + * + * cosh(z) = cosh(x+iy) + * = cosh(x) cos(y) + i sinh(x) sin(y). + * + * Exceptional values are noted in the comments within the source code. + * These values and the return value were taken from n1124.pdf. + * + * CCOSH_BIG is chosen such that + * spacing(0.5 * exp(CCOSH_BIG)) > 0.5*exp(-CCOSH_BIG) + * although the exact value assigned to CCOSH_BIG is not so important + */ +@ctype@ +npy_ccosh@c@(@ctype@ z) +{ +#if @precision@ == 1 + const npy_float CCOSH_BIG = 9.0f; + const npy_float CCOSH_HUGE = 1.70141183e+38f; +#endif +#if @precision@ == 2 + const npy_double CCOSH_BIG = 22.0; + const npy_double CCOSH_HUGE = 8.9884656743115795e+307; +#endif +#if @precision@ >= 3 +#if NPY_SIZEOF_LONGDOUBLE == NPY_SIZEOF_DOUBLE + const npy_longdouble CCOSH_BIG = 22.0L; + const npy_longdouble CCOSH_HUGE = 8.9884656743115795e+307L; +#else + const npy_longdouble CCOSH_BIG = 24.0L; + const npy_longdouble CCOSH_HUGE = 5.94865747678615882543e+4931L; +#endif +#endif + + @type@ x, y, h, absx; + npy_int xfinite, yfinite; + x = npy_creal@c@(z); y = npy_cimag@c@(z); - return npy_cpack@c@(npy_cos@c@(x) * npy_cosh@c@(y), -(npy_sin@c@(x) * npy_sinh@c@(y))); + + xfinite = npy_isfinite(x); + yfinite = npy_isfinite(y); + + /* Handle the nearly-non-exceptional cases where x and y are finite. */ + if (xfinite && yfinite) { + if (y == 0) { + return npy_cpack@c@(npy_cosh@c@(x), x * y); + } + absx = npy_fabs@c@(x); + if (absx < CCOSH_BIG) { + /* small x: normal case */ + return npy_cpack@c@(npy_cosh@c@(x) * npy_cos@c@(y), + npy_sinh@c@(x) * npy_sin@c@(y)); + } + + /* |x| >= 22, so cosh(x) ~= exp(|x|) */ + if (absx < SCALED_CEXP_LOWER@C@) { + /* x < 710: exp(|x|) won't overflow */ + h = npy_exp@c@(absx) * 0.5@c@; + return npy_cpack@c@(h * npy_cos@c@(y), + npy_copysign@c@(h, x) * npy_sin@c@(y)); + } + else if (absx < SCALED_CEXP_UPPER@C@) { + /* x < 1455: scale to avoid overflow */ + z = _npy_scaled_cexp@c@(absx, y, -1); + return npy_cpack@c@(npy_creal@c@(z), + npy_cimag@c@(z) * npy_copysign@c@(1, x)); + } + else { + /* x >= 1455: the result always overflows */ + h = CCOSH_HUGE * x; + return npy_cpack@c@(h * h * npy_cos@c@(y), h * npy_sin@c@(y)); + } + } + + /* + * cosh(+-0 +- I Inf) = dNaN + I sign(d(+-0, dNaN))0. + * The sign of 0 in the result is unspecified. Choice = normally + * the same as dNaN. Raise the invalid floating-point exception. + * + * cosh(+-0 +- I NaN) = d(NaN) + I sign(d(+-0, NaN))0. + * The sign of 0 in the result is unspecified. Choice = normally + * the same as d(NaN). + */ + if (x == 0 && !yfinite) { + return npy_cpack@c@(y - y, npy_copysign@c@(0, x * (y - y))); + } + + /* + * cosh(+-Inf +- I 0) = +Inf + I (+-)(+-)0. + * + * cosh(NaN +- I 0) = d(NaN) + I sign(d(NaN, +-0))0. + * The sign of 0 in the result is unspecified. + */ + if (y == 0 && !xfinite) { + return npy_cpack@c@(x * x, npy_copysign@c@(0, x) * y); + } + + /* + * cosh(x +- I Inf) = dNaN + I dNaN. + * Raise the invalid floating-point exception for finite nonzero x. + * + * cosh(x + I NaN) = d(NaN) + I d(NaN). + * Optionally raises the invalid floating-point exception for finite + * nonzero x. Choice = don't raise (except for signaling NaNs). + */ + if (xfinite && !yfinite) { + return npy_cpack@c@(y - y, x * (y - y)); + } + + /* + * cosh(+-Inf + I NaN) = +Inf + I d(NaN). + * + * cosh(+-Inf +- I Inf) = +Inf + I dNaN. + * The sign of Inf in the result is unspecified. Choice = always +. + * Raise the invalid floating-point exception. + * + * cosh(+-Inf + I y) = +Inf cos(y) +- I Inf sin(y) + */ + if (npy_isinf(x)) { + if (!yfinite) { + return npy_cpack@c@(x * x, x * (y - y)); + } + return npy_cpack@c@((x * x) * npy_cos@c@(y), x * npy_sin@c@(y)); + } + + /* + * cosh(NaN + I NaN) = d(NaN) + I d(NaN). + * + * cosh(NaN +- I Inf) = d(NaN) + I d(NaN). + * Optionally raises the invalid floating-point exception. + * Choice = raise. + * + * cosh(NaN + I y) = d(NaN) + I d(NaN). + * Optionally raises the invalid floating-point exception for finite + * nonzero y. Choice = don't raise (except for signaling NaNs). + */ + return npy_cpack@c@((x * x) * (y - y), (x + x) * (y - y)); } +#undef CCOSH_BIG +#undef CCOSH_HUGE #endif -#ifndef HAVE_CSIN@C@ -@ctype@ npy_csin@c@(@ctype@ z) +#ifndef HAVE_CSINH@C@ +/* + * Taken from the msun library in FreeBSD, rev 226599. + * + * Hyperbolic sine of a complex argument z = x + i y. + * + * sinh(z) = sinh(x+iy) + * = sinh(x) cos(y) + i cosh(x) sin(y). + * + * Exceptional values are noted in the comments within the source code. + * These values and the return value were taken from n1124.pdf. + */ +@ctype@ +npy_csinh@c@(@ctype@ z) +{ +#if @precision@ == 1 + const npy_float CSINH_BIG = 9.0f; + const npy_float CSINH_HUGE = 1.70141183e+38f; +#endif +#if @precision@ == 2 + const npy_double CSINH_BIG = 22.0; + const npy_double CSINH_HUGE = 8.9884656743115795e+307; +#endif +#if @precision@ >= 3 +#if NPY_SIZEOF_LONGDOUBLE == NPY_SIZEOF_DOUBLE + const npy_longdouble CSINH_BIG = 22.0L; + const npy_longdouble CSINH_HUGE = 8.9884656743115795e+307L; +#else + const npy_longdouble CSINH_BIG = 24.0L; + const npy_longdouble CSINH_HUGE = 5.94865747678615882543e+4931L; +#endif +#endif + + @type@ x, y, h, absx; + npy_int xfinite, yfinite; + + x = npy_creal@c@(z); + y = npy_cimag@c@(z); + + xfinite = npy_isfinite(x); + yfinite = npy_isfinite(y); + + /* Handle the nearly-non-exceptional cases where x and y are finite. */ + if (xfinite && yfinite) { + if (y == 0) { + return npy_cpack@c@(npy_sinh@c@(x), y); + } + absx = npy_fabs@c@(x); + if (absx < CSINH_BIG) { + /* small x: normal case */ + return npy_cpack@c@(npy_sinh@c@(x) * npy_cos@c@(y), + npy_cosh@c@(x) * npy_sin@c@(y)); + } + + /* |x| >= 22, so cosh(x) ~= exp(|x|) */ + if (absx < SCALED_CEXP_LOWER@C@) { + /* x < 710: exp(|x|) won't overflow */ + h = npy_exp@c@(npy_fabs@c@(x)) * 0.5@c@; + return npy_cpack@c@(npy_copysign@c@(h, x) * npy_cos@c@(y), + h * npy_sin@c@(y)); + } + else if (x < SCALED_CEXP_UPPER@C@) { + /* x < 1455: scale to avoid overflow */ + z = _npy_scaled_cexp@c@(absx, y, -1); + return npy_cpack@c@(npy_creal@c@(z) * npy_copysign@c@(1, x), + npy_cimag@c@(z)); + } + else { + /* x >= 1455: the result always overflows */ + h = CSINH_HUGE * x; + return npy_cpack@c@(h * npy_cos@c@(y), h * h * npy_sin@c@(y)); + } + } + + /* + * sinh(+-0 +- I Inf) = sign(d(+-0, dNaN))0 + I dNaN. + * The sign of 0 in the result is unspecified. Choice = normally + * the same as dNaN. Raise the invalid floating-point exception. + * + * sinh(+-0 +- I NaN) = sign(d(+-0, NaN))0 + I d(NaN). + * The sign of 0 in the result is unspecified. Choice = normally + * the same as d(NaN). + */ + if (x == 0 && !yfinite) { + return npy_cpack@c@(npy_copysign@c@(0, x * (y - y)), y - y); + } + + /* + * sinh(+-Inf +- I 0) = +-Inf + I +-0. + * + * sinh(NaN +- I 0) = d(NaN) + I +-0. + */ + if (y == 0 && !xfinite) { + if (npy_isnan(x)) { + return z; + } + return npy_cpack@c@(x, npy_copysign@c@(0, y)); + } + + /* + * sinh(x +- I Inf) = dNaN + I dNaN. + * Raise the invalid floating-point exception for finite nonzero x. + * + * sinh(x + I NaN) = d(NaN) + I d(NaN). + * Optionally raises the invalid floating-point exception for finite + * nonzero x. Choice = don't raise (except for signaling NaNs). + */ + if (xfinite && !yfinite) { + return npy_cpack@c@(y - y, x * (y - y)); + } + + /* + * sinh(+-Inf + I NaN) = +-Inf + I d(NaN). + * The sign of Inf in the result is unspecified. Choice = normally + * the same as d(NaN). + * + * sinh(+-Inf +- I Inf) = +Inf + I dNaN. + * The sign of Inf in the result is unspecified. Choice = always +. + * Raise the invalid floating-point exception. + * + * sinh(+-Inf + I y) = +-Inf cos(y) + I Inf sin(y) + */ + if (!xfinite && !npy_isnan(x)) { + if (!yfinite) { + return npy_cpack@c@(x * x, x * (y - y)); + } + return npy_cpack@c@(x * npy_cos@c@(y), + NPY_INFINITY@C@ * npy_sin@c@(y)); + } + + /* + * sinh(NaN + I NaN) = d(NaN) + I d(NaN). + * + * sinh(NaN +- I Inf) = d(NaN) + I d(NaN). + * Optionally raises the invalid floating-point exception. + * Choice = raise. + * + * sinh(NaN + I y) = d(NaN) + I d(NaN). + * Optionally raises the invalid floating-point exception for finite + * nonzero y. Choice = don't raise (except for signaling NaNs). + */ + return npy_cpack@c@((x * x) * (y - y), (x + x) * (y - y)); +} +#undef CSINH_BIG +#undef CSINH_HUGE +#endif + +#ifndef HAVE_CTANH@C@ +/* + * Taken from the msun library in FreeBSD, rev 226600. + * + * Hyperbolic tangent of a complex argument z = x + i y. + * + * The algorithm is from: + * + * W. Kahan. Branch Cuts for Complex Elementary Functions or Much + * Ado About Nothing's Sign Bit. In The State of the Art in + * Numerical Analysis, pp. 165 ff. Iserles and Powell, eds., 1987. + * + * Method: + * + * Let t = tan(x) + * beta = 1/cos^2(y) + * s = sinh(x) + * rho = cosh(x) + * + * We have: + * + * tanh(z) = sinh(z) / cosh(z) + * + * sinh(x) cos(y) + i cosh(x) sin(y) + * = --------------------------------- + * cosh(x) cos(y) + i sinh(x) sin(y) + * + * cosh(x) sinh(x) / cos^2(y) + i tan(y) + * = ------------------------------------- + * 1 + sinh^2(x) / cos^2(y) + * + * beta rho s + i t + * = ---------------- + * 1 + beta s^2 + * + * Modifications: + * + * I omitted the original algorithm's handling of overflow in tan(x) after + * verifying with nearpi.c that this can't happen in IEEE single or double + * precision. I also handle large x differently. + */ + +#define TANH_HUGE 22.0 +#define TANHF_HUGE 11.0F +#define TANHL_HUGE 42.0L + +@ctype@ +npy_ctanh@c@(@ctype@ z) { @type@ x, y; + @type@ t, beta, s, rho, denom; + + x = npy_creal@c@(z); + y = npy_cimag@c@(z); + + /* + * ctanh(NaN + i 0) = NaN + i 0 + * + * ctanh(NaN + i y) = NaN + i NaN for y != 0 + * + * The imaginary part has the sign of x*sin(2*y), but there's no + * special effort to get this right. + * + * ctanh(+-Inf +- i Inf) = +-1 +- 0 + * + * ctanh(+-Inf + i y) = +-1 + 0 sin(2y) for y finite + * + * The imaginary part of the sign is unspecified. This special + * case is only needed to avoid a spurious invalid exception when + * y is infinite. + */ + if (!npy_isfinite(x)) { + if (npy_isnan(x)) { + return npy_cpack@c@(x, (y == 0 ? y : x * y)); + } + return npy_cpack@c@(npy_copysign@c@(1,x), + npy_copysign@c@(0, + npy_isinf(y) ? + y : npy_sin@c@(y) * npy_cos@c@(y))); + } + + /* + * ctanh(x + i NAN) = NaN + i NaN + * ctanh(x +- i Inf) = NaN + i NaN + */ + if (!npy_isfinite(y)) { + return (npy_cpack@c@(y - y, y - y)); + } + + /* + * ctanh(+-huge + i +-y) ~= +-1 +- i 2sin(2y)/exp(2x), using the + * approximation sinh^2(huge) ~= exp(2*huge) / 4. + * We use a modified formula to avoid spurious overflow. + */ + if (npy_fabs@c@(x) >= TANH@C@_HUGE) { + @type@ exp_mx = npy_exp@c@(-npy_fabs@c@(x)); + return npy_cpack@c@(npy_copysign@c@(1, x), + 4 * npy_sin@c@(y) * npy_cos@c@(y) * + exp_mx * exp_mx); + } + + /* Kahan's algorithm */ + t = npy_tan@c@(y); + beta = 1 + t * t; /* = 1 / cos^2(y) */ + s = npy_sinh@c@(x); + rho = npy_sqrt@c@(1 + s * s); /* = cosh(x) */ + denom = 1 + beta * s * s; + return (npy_cpack@c@((beta * rho * s) / denom, t / denom)); +} +#undef TANH_HUGE +#undef TANHF_HUGE +#undef TANHL_HUGE +#endif + +#if !defined (HAVE_CACOS@C@) || !defined (HAVE_CASINH@C@) +/* + * Complex inverse trig functions taken from the msum library in FreeBSD + * revision 251404 + * + * The algorithm is very close to that in "Implementing the complex arcsine + * and arccosine functions using exception handling" by T. E. Hull, Thomas F. + * Fairgrieve, and Ping Tak Peter Tang, published in ACM Transactions on + * Mathematical Software, Volume 23 Issue 3, 1997, Pages 299-335, + * http://dl.acm.org/citation.cfm?id=275324. + * + * Throughout we use the convention z = x + I*y. + * + * casinh(z) = sign(x)*log(A+sqrt(A*A-1)) + I*asin(B) + * where + * A = (|z+I| + |z-I|) / 2 + * B = (|z+I| - |z-I|) / 2 = y/A + * + * These formulas become numerically unstable: + * (a) for Re(casinh(z)) when z is close to the line segment [-I, I] (that + * is, Re(casinh(z)) is close to 0); + * (b) for Im(casinh(z)) when z is close to either of the intervals + * [I, I*infinity) or (-I*infinity, -I] (that is, |Im(casinh(z))| is + * close to PI/2). + * + * These numerical problems are overcome by defining + * f(a, b) = (hypot(a, b) - b) / 2 = a*a / (hypot(a, b) + b) / 2 + * Then if A < A_crossover, we use + * log(A + sqrt(A*A-1)) = log1p((A-1) + sqrt((A-1)*(A+1))) + * A-1 = f(x, 1+y) + f(x, 1-y) + * and if B > B_crossover, we use + * asin(B) = atan2(y, sqrt(A*A - y*y)) = atan2(y, sqrt((A+y)*(A-y))) + * A-y = f(x, y+1) + f(x, y-1) + * where without loss of generality we have assumed that x and y are + * non-negative. + * + * Much of the difficulty comes because the intermediate computations may + * produce overflows or underflows. This is dealt with in the paper by Hull + * et al by using exception handling. We do this by detecting when + * computations risk underflow or overflow. The hardest part is handling the + * underflows when computing f(a, b). + * + * Note that the function f(a, b) does not appear explicitly in the paper by + * Hull et al, but the idea may be found on pages 308 and 309. Introducing the + * function f(a, b) allows us to concentrate many of the clever tricks in this + * paper into one function. + */ + +/* + * Function f(a, b, hypot_a_b) = (hypot(a, b) - b) / 2. + * Pass hypot(a, b) as the third argument. + */ +static NPY_INLINE @type@ +_f@c@(@type@ a, @type@ b, @type@ hypot_a_b) +{ + if (b < 0) { + return ((hypot_a_b - b) / 2); + } + if (b == 0) { + return (a / 2); + } + return (a * a / (hypot_a_b + b) / 2); +} + +/* + * All the hard work is contained in this function. + * x and y are assumed positive or zero, and less than RECIP_EPSILON. + * Upon return: + * rx = Re(casinh(z)) = -Im(cacos(y + I*x)). + * B_is_usable is set to 1 if the value of B is usable. + * If B_is_usable is set to 0, sqrt_A2my2 = sqrt(A*A - y*y), and new_y = y. + * If returning sqrt_A2my2 has potential to result in an underflow, it is + * rescaled, and new_y is similarly rescaled. + */ +static NPY_INLINE void +_do_hard_work@c@(@type@ x, @type@ y, @type@ *rx, + npy_int *B_is_usable, @type@ *B, @type@ *sqrt_A2my2, @type@ *new_y) +{ +#if @precision@ == 1 + const npy_float A_crossover = 10.0f; + const npy_float B_crossover = 0.6417f; + const npy_float FOUR_SQRT_MIN = 4.3368086899420177e-19f; +#endif +#if @precision@ == 2 + const npy_double A_crossover = 10.0; + const npy_double B_crossover = 0.6417; + const npy_double FOUR_SQRT_MIN = 5.9666725849601654e-154; +#endif +#if @precision@ == 3 + const npy_longdouble A_crossover = 10.0l; + const npy_longdouble B_crossover = 0.6417l; +#if NPy_SIZEOF_LONGDOUBLE == NPY_SIZEOF_DOUBLE + const npy_longdouble FOUR_SQRT_MIN = 5.9666725849601654e-154; +#else + const npy_longdouble FOUR_SQRT_MIN = 7.3344154702193886625e-2466l; +#endif +#endif + @type@ R, S, A; /* A, B, R, and S are as in Hull et al. */ + @type@ Am1, Amy; /* A-1, A-y. */ + + R = npy_hypot@c@(x, y + 1); /* |z+I| */ + S = npy_hypot@c@(x, y - 1); /* |z-I| */ + + /* A = (|z+I| + |z-I|) / 2 */ + A = (R + S) / 2; + /* + * Mathematically A >= 1. There is a small chance that this will not + * be so because of rounding errors. So we will make certain it is + * so. + */ + if (A < 1) { + A = 1; + } + + if (A < A_crossover) { + /* + * Am1 = fp + fm, where fp = f(x, 1+y), and fm = f(x, 1-y). + * rx = log1p(Am1 + sqrt(Am1*(A+1))) + */ + if (y == 1 && x < @TEPS@ * @TEPS@ / 128) { + /* + * fp is of order x^2, and fm = x/2. + * A = 1 (inexactly). + */ + *rx = npy_sqrt@c@(x); + } + else if (x >= @TEPS@ * npy_fabs@c@(y - 1)) { + /* + * Underflow will not occur because + * x >= DBL_EPSILON^2/128 >= FOUR_SQRT_MIN + */ + Am1 = _f@c@(x, 1 + y, R) + _f@c@(x, 1 - y, S); + *rx = npy_log1p@c@(Am1 + npy_sqrt@c@(Am1 * (A + 1))); + } + else if (y < 1) { + /* + * fp = x*x/(1+y)/4, fm = x*x/(1-y)/4, and + * A = 1 (inexactly). + */ + *rx = x / npy_sqrt@c@((1 - y) * (1 + y)); + } + else { /* if (y > 1) */ + /* + * A-1 = y-1 (inexactly). + */ + *rx = npy_log1p@c@((y - 1) + npy_sqrt@c@((y - 1) * (y + 1))); + } + } + else { + *rx = npy_log@c@(A + npy_sqrt@c@(A * A - 1)); + } + + *new_y = y; + + if (y < FOUR_SQRT_MIN) { + /* + * Avoid a possible underflow caused by y/A. For casinh this + * would be legitimate, but will be picked up by invoking atan2 + * later on. For cacos this would not be legitimate. + */ + *B_is_usable = 0; + *sqrt_A2my2 = A * (2 / @TEPS@); + *new_y = y * (2 / @TEPS@); + return; + } + + /* B = (|z+I| - |z-I|) / 2 = y/A */ + *B = y / A; + *B_is_usable = 1; + + if (*B > B_crossover) { + *B_is_usable = 0; + /* + * Amy = fp + fm, where fp = f(x, y+1), and fm = f(x, y-1). + * sqrt_A2my2 = sqrt(Amy*(A+y)) + */ + if (y == 1 && x < @TEPS@ / 128) { + /* + * fp is of order x^2, and fm = x/2. + * A = 1 (inexactly). + */ + *sqrt_A2my2 = npy_sqrt@c@(x) * npy_sqrt@c@((A + y) / 2); + } + else if (x >= @TEPS@ * npy_fabs@c@(y - 1)) { + /* + * Underflow will not occur because + * x >= DBL_EPSILON/128 >= FOUR_SQRT_MIN + * and + * x >= DBL_EPSILON^2 >= FOUR_SQRT_MIN + */ + Amy = _f@c@(x, y + 1, R) + _f@c@(x, y - 1, S); + *sqrt_A2my2 = npy_sqrt@c@(Amy * (A + y)); + } + else if (y > 1) { + /* + * fp = x*x/(y+1)/4, fm = x*x/(y-1)/4, and + * A = y (inexactly). + * + * y < RECIP_EPSILON. So the following + * scaling should avoid any underflow problems. + */ + *sqrt_A2my2 = x * (4 / @TEPS@ / @TEPS@) * y / + npy_sqrt@c@((y + 1) * (y - 1)); + *new_y = y * (4 / @TEPS@ / @TEPS@); + } + else { /* if (y < 1) */ + /* + * fm = 1-y >= DBL_EPSILON, fp is of order x^2, and + * A = 1 (inexactly). + */ + *sqrt_A2my2 = npy_sqrt@c@((1 - y) * (1 + y)); + } + } +} + +/* + * Optimized version of clog() for |z| finite and larger than ~RECIP_EPSILON. + */ +static NPY_INLINE void +_clog_for_large_values@c@(@type@ x, @type@ y, + @type@ *rr, @type@ *ri) +{ +#if @precision@ == 1 + const npy_float QUARTER_SQRT_MAX = 4.611685743549481e+18f; + const npy_float SQRT_MIN = 1.0842021724855044e-19f; + #endif +#if @precision@ == 2 + const npy_double QUARTER_SQRT_MAX = 3.3519519824856489e+153; + const npy_double SQRT_MIN = 1.4916681462400413e-154; + #endif +#if @precision@ == 3 +#if NPY_SIZEOF_LONGDOUBLE == NPY_SIZEOF_DOUBLE + const npy_longdouble QUARTER_SQRT_MAX = 3.3519519824856489e+153; + const npy_longdouble SQRT_MIN = 1.4916681462400413e-154; +#else + const npy_longdouble QUARTER_SQRT_MAX = 2.7268703390485398235e+2465l; + const npy_longdouble SQRT_MIN = 1.8336038675548471656e-2466l; +#endif +#endif + @type@ ax, ay, t; + + ax = npy_fabs@c@(x); + ay = npy_fabs@c@(y); + if (ax < ay) { + t = ax; + ax = ay; + ay = t; + } + + /* + * Avoid overflow in hypot() when x and y are both very large. + * Divide x and y by E, and then add 1 to the logarithm. This depends + * on E being larger than sqrt(2). + * Dividing by E causes an insignificant loss of accuracy; however + * this method is still poor since it is uneccessarily slow. + */ + if (ax > @TMAX@ / 2) { + *rr = npy_log@c@(npy_hypot@c@(x / NPY_E@c@, y / NPY_E@c@)) + 1; + } + /* + * Avoid overflow when x or y is large. Avoid underflow when x or + * y is small. + */ + else if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN) { + *rr = npy_log@c@(npy_hypot@c@(x, y)); + } + else { + *rr = npy_log@c@(ax * ax + ay * ay) / 2; + } + *ri = npy_atan2@c@(y, x); +} +#endif + +#ifndef HAVE_CACOS@C@ +@ctype@ +npy_cacos@c@(@ctype@ z) +{ +#if @precision@ == 1 + /* this is sqrt(6*EPS) */ + const npy_float SQRT_6_EPSILON = 8.4572793338e-4f; + /* chosen such that pio2_hi + pio2_lo == pio2_hi but causes FE_INEXACT. */ + const volatile npy_float pio2_lo = 7.5497899549e-9f; +#endif +#if @precision@ == 2 + const npy_double SQRT_6_EPSILON = 3.65002414998885671e-08; + const volatile npy_double pio2_lo = 6.1232339957367659e-17; +#endif +#if @precision@ == 3 + const npy_longdouble SQRT_6_EPSILON = 8.0654900873493277169e-10l; + const volatile npy_longdouble pio2_lo = 2.710505431213761085e-20l; +#endif + const @type@ RECIP_EPSILON = 1.0@c@ / @TEPS@; + const @type@ pio2_hi = NPY_PI_2@c@; + @type@ x, y, ax, ay, wx, wy, rx, ry, B, sqrt_A2mx2, new_x; + npy_int sx, sy; + npy_int B_is_usable; + x = npy_creal@c@(z); y = npy_cimag@c@(z); - return npy_cpack@c@(npy_sin@c@(x) * npy_cosh@c@(y), npy_cos@c@(x) * npy_sinh@c@(y)); + sx = npy_signbit(x); + sy = npy_signbit(y); + ax = npy_fabs@c@(x); + ay = npy_fabs@c@(y); + + if (npy_isnan(x) || npy_isnan(y)) { + /* cacos(+-Inf + I*NaN) = NaN + I*opt(-)Inf */ + if (npy_isinf(x)) { + return npy_cpack@c@(y + y, -NPY_INFINITY@C@); + } + /* cacos(NaN + I*+-Inf) = NaN + I*-+Inf */ + if (npy_isinf(y)) { + return npy_cpack@c@(x + x, -y); + } + /* cacos(0 + I*NaN) = PI/2 + I*NaN with inexact */ + if (x == 0) { + return npy_cpack@c@(pio2_hi + pio2_lo, y + y); + } + /* + * All other cases involving NaN return NaN + I*NaN. + * C99 leaves it optional whether to raise invalid if one of + * the arguments is not NaN, so we opt not to raise it. + */ + return npy_cpack@c@(NPY_NAN@C@, NPY_NAN@C@); + } + + if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { + /* clog...() will raise inexact unless x or y is infinite. */ + _clog_for_large_values@c@(x, y, &wx, &wy); + rx = npy_fabs@c@(wy); + ry = wx + NPY_LOGE2@c@; + if (sy == 0) { + ry = -ry; + } + return npy_cpack@c@(rx, ry); + } + + /* Avoid spuriously raising inexact for z = 1. */ + if (x == 1 && y == 0) { + return npy_cpack@c@(0, -y); + } + + /* All remaining cases are inexact. */ + raise_inexact(); + + if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4) { + return npy_cpack@c@(pio2_hi - (x - pio2_lo), -y); + } + + _do_hard_work@c@(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x); + if (B_is_usable) { + if (sx == 0) { + rx = npy_acos@c@(B); + } + else { + rx = npy_acos@c@(-B); + } + } + else { + if (sx == 0) { + rx = npy_atan2@c@(sqrt_A2mx2, new_x); + } + else { + rx = npy_atan2@c@(sqrt_A2mx2, -new_x); + } + } + if (sy == 0) { + ry = -ry; + } + return npy_cpack@c@(rx, ry); +} +#endif + +#ifndef HAVE_CASIN@C@ +@ctype@ +npy_casin@c@(@ctype@ z) +{ + /* casin(z) = I * conj( casinh(I * conj(z)) ) */ + z = npy_casinh@c@(npy_cpack@c@(npy_cimag@c@(z), npy_creal@c@(z))); + return npy_cpack@c@(npy_cimag@c@(z), npy_creal@c@(z)); +} +#endif + +#ifndef HAVE_CATAN@C@ +@ctype@ +npy_catan@c@(@ctype@ z) +{ + /* catan(z) = I * conj( catanh(I * conj(z)) ) */ + z = npy_catanh@c@(npy_cpack@c@(npy_cimag@c@(z), npy_creal@c@(z))); + return npy_cpack@c@(npy_cimag@c@(z), npy_creal@c@(z)); +} +#endif + +#ifndef HAVE_CACOSH@C@ +@ctype@ +npy_cacosh@c@(@ctype@ z) +{ + /* + * cacosh(z) = I*cacos(z) or -I*cacos(z) + * where the sign is chosen so Re(cacosh(z)) >= 0. + */ + @ctype@ w; + @type@ rx, ry; + + w = npy_cacos@c@(z); + rx = npy_creal@c@(w); + ry = npy_cimag@c@(w); + /* cacosh(NaN + I*NaN) = NaN + I*NaN */ + if (npy_isnan(rx) && npy_isnan(ry)) { + return npy_cpack@c@(ry, rx); + } + /* cacosh(NaN + I*+-Inf) = +Inf + I*NaN */ + /* cacosh(+-Inf + I*NaN) = +Inf + I*NaN */ + if (npy_isnan(rx)) { + return npy_cpack@c@(npy_fabs@c@(ry), rx); + } + /* cacosh(0 + I*NaN) = NaN + I*NaN */ + if (npy_isnan(ry)) { + return npy_cpack@c@(ry, ry); + } + return npy_cpack@c@(npy_fabs@c@(ry), npy_copysign@c@(rx, npy_cimag@c@(z))); +} +#endif + +#ifndef HAVE_CASINH@C@ +/* + * casinh(z) = z + O(z^3) as z -> 0 + * + * casinh(z) = sign(x)*clog(sign(x)*z) + O(1/z^2) as z -> infinity + * The above formula works for the imaginary part as well, because + * Im(casinh(z)) = sign(x)*atan2(sign(x)*y, fabs(x)) + O(y/z^3) + * as z -> infinity, uniformly in y + */ +@ctype@ +npy_casinh@c@(@ctype@ z) +{ +#if @precision@ == 1 + /* this is sqrt(6*EPS) */ + const npy_float SQRT_6_EPSILON = 8.4572793338e-4f; + /* chosen such that pio2_hi + pio2_lo == pio2_hi but causes FE_INEXACT. */ + const volatile npy_float pio2_lo = 7.5497899549e-9f; +#endif +#if @precision@ == 2 + const npy_double SQRT_6_EPSILON = 3.65002414998885671e-08; + const volatile npy_double pio2_lo = 6.1232339957367659e-17; +#endif +#if @precision@ == 3 + const npy_longdouble SQRT_6_EPSILON = 8.0654900873493277169e-10l; + const volatile npy_longdouble pio2_lo = 2.710505431213761085e-20l; +#endif + const @type@ RECIP_EPSILON = 1.0@c@ / @TEPS@; + const @type@ pio2_hi = NPY_PI_2@c@; + @type@ x, y, ax, ay, wx, wy, rx, ry, B, sqrt_A2my2, new_y; + npy_int B_is_usable; + + x = npy_creal@c@(z); + y = npy_cimag@c@(z); + ax = npy_fabs@c@(x); + ay = npy_fabs@c@(y); + + if (npy_isnan(x) || npy_isnan(y)) { + /* casinh(+-Inf + I*NaN) = +-Inf + I*NaN */ + if (npy_isinf(x)) { + return npy_cpack@c@(x, y + y); + } + /* casinh(NaN + I*+-Inf) = opt(+-)Inf + I*NaN */ + if (npy_isinf(y)) { + return npy_cpack@c@(y, x + x); + } + /* casinh(NaN + I*0) = NaN + I*0 */ + if (y == 0) { + return npy_cpack@c@(x + x, y); + } + /* + * All other cases involving NaN return NaN + I*NaN. + * C99 leaves it optional whether to raise invalid if one of + * the arguments is not NaN, so we opt not to raise it. + */ + return npy_cpack@c@(NPY_NAN@C@, NPY_NAN@C@); + } + + if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { + /* clog...() will raise inexact unless x or y is infinite. */ + if (npy_signbit(x) == 0) { + _clog_for_large_values@c@(x, y, &wx, &wy); + wx += NPY_LOGE2@c@; + } + else { + _clog_for_large_values@c@(-x, -y, &wx, &wy); + wx += NPY_LOGE2@c@; + } + return npy_cpack@c@(npy_copysign@c@(wx, x), npy_copysign@c@(wy, y)); + } + + /* Avoid spuriously raising inexact for z = 0. */ + if (x == 0 && y == 0) { + return (z); + } + + /* All remaining cases are inexact. */ + raise_inexact(); + + if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4) { + return (z); + } + + _do_hard_work@c@(ax, ay, &rx, &B_is_usable, &B, &sqrt_A2my2, &new_y); + if (B_is_usable) { + ry = npy_asin@c@(B); + } + else { + ry = npy_atan2@c@(new_y, sqrt_A2my2); + } + return npy_cpack@c@(npy_copysign@c@(rx, x), npy_copysign@c@(ry, y)); +} +#endif + +#ifndef HAVE_CATANH@C@ +/* + * sum_squares(x,y) = x*x + y*y (or just x*x if y*y would underflow). + * Assumes x*x and y*y will not overflow. + * Assumes x and y are finite. + * Assumes y is non-negative. + * Assumes fabs(x) >= DBL_EPSILON. + */ +static NPY_INLINE @type@ +_sum_squares@c@(@type@ x, @type@ y) +{ +#if @precision@ == 1 +const npy_float SQRT_MIN = 1.0842022e-19f; +#endif +#if @precision@ == 2 +const npy_double SQRT_MIN = 1.4916681462400413e-154; /* sqrt(DBL_MIN) */ +#endif +#if @precision@ == 3 +/* this is correct for 80 bit long doubles */ +const npy_longdouble SQRT_MIN = 1.8336038675548471656e-2466l; +#endif + /* Avoid underflow when y is small. */ + if (y < SQRT_MIN) { + return (x * x); + } + + return (x * x + y * y); +} + +/* + * real_part_reciprocal(x, y) = Re(1/(x+I*y)) = x/(x*x + y*y). + * Assumes x and y are not NaN, and one of x and y is larger than + * RECIP_EPSILON. We avoid unwarranted underflow. It is important to not use + * the code creal(1/z), because the imaginary part may produce an unwanted + * underflow. + * This is only called in a context where inexact is always raised before + * the call, so no effort is made to avoid or force inexact. + */ +#if @precision@ == 1 +#define BIAS (FLT_MAX_EXP - 1) +#define CUTOFF (FLT_MANT_DIG / 2 + 1) +static NPY_INLINE npy_float +_real_part_reciprocalf(npy_float x, npy_float y) +{ + npy_float scale; + npy_uint32 hx, hy; + npy_int32 ix, iy; + + GET_FLOAT_WORD(hx, x); + ix = hx & 0x7f800000; + GET_FLOAT_WORD(hy, y); + iy = hy & 0x7f800000; + if (ix - iy >= CUTOFF << 23 || npy_isinf(x)) { + return (1 / x); + } + if (iy - ix >= CUTOFF << 23) { + return (x / y / y); + } + if (ix <= (BIAS + FLT_MAX_EXP / 2 - CUTOFF) << 23) { + return (x / (x * x + y * y)); + } + SET_FLOAT_WORD(scale, 0x7f800000 - ix); + x *= scale; + y *= scale; + return (x / (x * x + y * y) * scale); +} +#undef BIAS +#undef CUTOFF +#endif +#if @precision@ == 2 +#define BIAS (DBL_MAX_EXP - 1) +/* more guard digits are useful iff there is extra precision. */ +#define CUTOFF (DBL_MANT_DIG / 2 + 1) /* just half or 1 guard digit */ +static NPY_INLINE npy_double +_real_part_reciprocal(npy_double x, npy_double y) +{ + npy_double scale; + npy_uint32 hx, hy; + npy_int32 ix, iy; + + /* + * This code is inspired by the C99 document n1124.pdf, Section G.5.1, + * example 2. + */ + GET_HIGH_WORD(hx, x); + ix = hx & 0x7ff00000; + GET_HIGH_WORD(hy, y); + iy = hy & 0x7ff00000; + if (ix - iy >= CUTOFF << 20 || npy_isinf(x)) { + /* +-Inf -> +-0 is special */ + return (1 / x); + } + if (iy - ix >= CUTOFF << 20) { + /* should avoid double div, but hard */ + return (x / y / y); + } + if (ix <= (BIAS + DBL_MAX_EXP / 2 - CUTOFF) << 20) { + return (x / (x * x + y * y)); + } + scale = 1; + SET_HIGH_WORD(scale, 0x7ff00000 - ix); /* 2**(1-ilogb(x)) */ + x *= scale; + y *= scale; + return (x / (x * x + y * y) * scale); +} +#undef BIAS +#undef CUTOFF +#endif +#if @precision@ == 3 +#ifndef HAVE_LDOUBLE_DOUBLE_DOUBLE_BE +#define BIAS (LDBL_MAX_EXP - 1) +#define CUTOFF (LDBL_MANT_DIG / 2 + 1) +static NPY_INLINE npy_longdouble +_real_part_reciprocall(npy_longdouble x, + npy_longdouble y) +{ + npy_longdouble scale; + union IEEEl2bitsrep ux, uy, us; + npy_int32 ix, iy; + + ux.e = x; + ix = GET_LDOUBLE_EXP(ux); + uy.e = y; + iy = GET_LDOUBLE_EXP(uy); + if (ix - iy >= CUTOFF || npy_isinf(x)) { + return (1/x); + } + if (iy - ix >= CUTOFF) { + return (x/y/y); + } + if (ix <= BIAS + LDBL_MAX_EXP / 2 - CUTOFF) { + return (x/(x*x + y*y)); + } + us.e = 1; + SET_LDOUBLE_EXP(us, 0x7fff - ix); + scale = us.e; + x *= scale; + y *= scale; + return (x/(x*x + y*y) * scale); +} +#undef BIAS +#undef CUTOFF +#else +static NPY_INLINE npy_longdouble +_real_part_reciprocall(npy_longdouble x, + npy_longdouble y) +{ + return x/(x*x + y*y); +} +#endif +#endif + +@ctype@ +npy_catanh@c@(@ctype@ z) +{ +#if @precision@ == 1 + /* this is sqrt(3*EPS) */ + const npy_float SQRT_3_EPSILON = 5.9801995673e-4f; + /* chosen such that pio2_hi + pio2_lo == pio2_hi but causes FE_INEXACT. */ + const volatile npy_float pio2_lo = 7.5497899549e-9f; +#endif +#if @precision@ == 2 + const npy_double SQRT_3_EPSILON = 2.5809568279517849e-8; + const volatile npy_double pio2_lo = 6.1232339957367659e-17; +#endif +#if @precision@ == 3 + const npy_longdouble SQRT_3_EPSILON = 5.70316273435758915310e-10l; + const volatile npy_longdouble pio2_lo = 2.710505431213761085e-20l; +#endif + const @type@ RECIP_EPSILON = 1.0@c@ / @TEPS@; + const @type@ pio2_hi = NPY_PI_2@c@; + @type@ x, y, ax, ay, rx, ry; + + x = npy_creal@c@(z); + y = npy_cimag@c@(z); + ax = npy_fabs@c@(x); + ay = npy_fabs@c@(y); + + /* This helps handle many cases. */ + if (y == 0 && ax <= 1) { + return npy_cpack@c@(npy_atanh@c@(x), y); + } + + /* To ensure the same accuracy as atan(), and to filter out z = 0. */ + if (x == 0) { + return npy_cpack@c@(x, npy_atan@c@(y)); + } + + if (npy_isnan(x) || npy_isnan(y)) { + /* catanh(+-Inf + I*NaN) = +-0 + I*NaN */ + if (npy_isinf(x)) { + return npy_cpack@c@(npy_copysign@c@(0, x), y + y); + } + /* catanh(NaN + I*+-Inf) = sign(NaN)0 + I*+-PI/2 */ + if (npy_isinf(y)) { + return npy_cpack@c@(npy_copysign@c@(0, x), + npy_copysign@c@(pio2_hi + pio2_lo, y)); + } + /* + * All other cases involving NaN return NaN + I*NaN. + * C99 leaves it optional whether to raise invalid if one of + * the arguments is not NaN, so we opt not to raise it. + */ + return npy_cpack@c@(NPY_NAN@C@, NPY_NAN@C@); + } + + if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { + return npy_cpack@c@(_real_part_reciprocal@c@(x, y), + npy_copysign@c@(pio2_hi + pio2_lo, y)); + } + + if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) { + /* + * z = 0 was filtered out above. All other cases must raise + * inexact, but this is the only only that needs to do it + * explicitly. + */ + raise_inexact(); + return (z); + } + + if (ax == 1 && ay < @TEPS@) { + rx = (NPY_LOGE2@c@ - npy_log@c@(ay)) / 2; + } + else { + rx = npy_log1p@c@(4 * ax / _sum_squares@c@(ax - 1, ay)) / 4; + } + + if (ax == 1) { + ry = npy_atan2@c@(2, -ay) / 2; + } + else if (ay < @TEPS@) { + ry = npy_atan2@c@(2 * ay, (1 - ax) * (1 + ax)) / 2; + } + else { + ry = npy_atan2@c@(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2; + } + + return npy_cpack@c@(npy_copysign@c@(rx, x), npy_copysign@c@(ry, y)); } #endif /**end repeat**/ @@ -245,7 +1753,8 @@ * #KIND = CABS,CARG# */ #ifdef HAVE_@KIND@@C@ -@type@ npy_@kind@@c@(@ctype@ z) +@type@ +npy_@kind@@c@(@ctype@ z) { __@ctype@_to_c99_cast z1; z1.npy_z = z; @@ -255,11 +1764,14 @@ /**end repeat1**/ /**begin repeat1 - * #kind = cexp,clog,csqrt,ccos,csin# - * #KIND = CEXP,CLOG,CSQRT,CCOS,CSIN# + * #kind = cexp,clog,csqrt,ccos,csin,ctan,ccosh,csinh,ctanh, + * cacos,casin,catan,cacosh,casinh,catanh# + * #KIND = CEXP,CLOG,CSQRT,CCOS,CSIN,CTAN,CCOSH,CSINH,CTANH, + * CACOS,CASIN,CATAN,CACOSH,CASINH,CATANH# */ #ifdef HAVE_@KIND@@C@ -@ctype@ npy_@kind@@c@(@ctype@ z) +@ctype@ +npy_@kind@@c@(@ctype@ z) { __@ctype@_to_c99_cast z1; __@ctype@_to_c99_cast ret; @@ -270,22 +1782,5 @@ #endif /**end repeat1**/ -/**begin repeat1 - * #kind = cpow# - * #KIND = CPOW# - */ -#ifdef HAVE_@KIND@@C@ -@ctype@ npy_@kind@@c@(@ctype@ x, @ctype@ y) -{ - __@ctype@_to_c99_cast xcast; - __@ctype@_to_c99_cast ycast; - __@ctype@_to_c99_cast ret; - xcast.npy_z = x; - ycast.npy_z = y; - ret.c99_z = @kind@@c@(xcast.c99_z, ycast.c99_z); - return ret.npy_z; -} -#endif -/**end repeat1**/ /**end repeat**/ diff --git a/numpy/core/src/private/npy_config.h b/numpy/core/src/private/npy_config.h index 882913e2f..6e98dc7e9 100644 --- a/numpy/core/src/private/npy_config.h +++ b/numpy/core/src/private/npy_config.h @@ -4,12 +4,6 @@ #include "config.h" #include "numpy/numpyconfig.h" -/* Disable broken MS math functions */ -#if defined(_MSC_VER) || defined(__MINGW32_VERSION) -#undef HAVE_ATAN2 -#undef HAVE_HYPOT -#endif - /* * largest alignment the copy loops might require * required as string, void and complex types might get copied using larger @@ -21,9 +15,56 @@ */ #define NPY_MAX_COPY_ALIGNMENT 16 +/* blacklist */ + /* Disable broken Sun Workshop Pro math functions */ #ifdef __SUNPRO_C + +#undef HAVE_ATAN2 +#undef HAVE_ATAN2F +#undef HAVE_ATAN2L + +#endif + +/* Disable broken MS math functions */ +#if defined(_MSC_VER) || defined(__MINGW32_VERSION) + #undef HAVE_ATAN2 +#undef HAVE_ATAN2F +#undef HAVE_ATAN2L + +#undef HAVE_HYPOT +#undef HAVE_HYPOTF +#undef HAVE_HYPOTL + #endif +/* Disable broken gnu trig functions on linux */ +#if defined(__linux__) && defined(__GNUC__) + +#if defined(HAVE_FEATURES_H) +#include <features.h> +#define TRIG_OK __GLIBC_PREREQ(2, 16) +#else +#define TRIG_OK 0 +#endif + +#if !TRIG_OK +#undef HAVE_CASIN +#undef HAVE_CASINF +#undef HAVE_CASINL +#undef HAVE_CASINH +#undef HAVE_CASINHF +#undef HAVE_CASINHL +#undef HAVE_CATAN +#undef HAVE_CATANF +#undef HAVE_CATANL +#undef HAVE_CATANH +#undef HAVE_CATANHF +#undef HAVE_CATANHL +#endif +#undef TRIG_OK + +#endif /* defined(__linux__) && defined(__GNUC__) */ + #endif diff --git a/numpy/core/src/umath/funcs.inc.src b/numpy/core/src/umath/funcs.inc.src index 3aad44c9f..f0fa3e9dd 100644 --- a/numpy/core/src/umath/funcs.inc.src +++ b/numpy/core/src/umath/funcs.inc.src @@ -177,49 +177,8 @@ npy_ObjectLogicalNot(PyObject *i1) * #ctype = npy_cfloat, npy_cdouble, npy_clongdouble# * #ftype = npy_float, npy_double, npy_longdouble# * #c = f, ,l# - * #C = F, ,L# - * #precision = 1,2,4# */ -/* - * Perform the operation result := 1 + coef * x * result, - * with real coefficient `coef`. - */ -#define SERIES_HORNER_TERM@c@(result, x, coef) \ - do { \ - nc_prod@c@((result), (x), (result)); \ - (result)->real *= (coef); \ - (result)->imag *= (coef); \ - nc_sum@c@((result), &nc_1@c@, (result)); \ - } while(0) - -/* constants */ -static @ctype@ nc_1@c@ = {1., 0.}; -static @ctype@ nc_half@c@ = {0.5, 0.}; -static @ctype@ nc_i@c@ = {0., 1.}; -static @ctype@ nc_i2@c@ = {0., 0.5}; -/* - * static @ctype@ nc_mi@c@ = {0.0@c@, -1.0@c@}; - * static @ctype@ nc_pi2@c@ = {NPY_PI_2@c@., 0.0@c@}; - */ - - -static void -nc_sum@c@(@ctype@ *a, @ctype@ *b, @ctype@ *r) -{ - r->real = a->real + b->real; - r->imag = a->imag + b->imag; - return; -} - -static void -nc_diff@c@(@ctype@ *a, @ctype@ *b, @ctype@ *r) -{ - r->real = a->real - b->real; - r->imag = a->imag - b->imag; - return; -} - static void nc_neg@c@(@ctype@ *a, @ctype@ *r) { @@ -229,26 +188,6 @@ nc_neg@c@(@ctype@ *a, @ctype@ *r) } static void -nc_prod@c@(@ctype@ *a, @ctype@ *b, @ctype@ *r) -{ - @ftype@ ar=a->real, br=b->real, ai=a->imag, bi=b->imag; - r->real = ar*br - ai*bi; - r->imag = ar*bi + ai*br; - return; -} - -static void -nc_quot@c@(@ctype@ *a, @ctype@ *b, @ctype@ *r) -{ - - @ftype@ ar=a->real, br=b->real, ai=a->imag, bi=b->imag; - @ftype@ d = br*br + bi*bi; - r->real = (ar*br + ai*bi)/d; - r->imag = (ai*br - ar*bi)/d; - return; -} - -static void nc_sqrt@c@(@ctype@ *x, @ctype@ *r) { *r = npy_csqrt@c@(*x); @@ -307,164 +246,28 @@ nc_expm1@c@(@ctype@ *x, @ctype@ *r) static void nc_pow@c@(@ctype@ *a, @ctype@ *b, @ctype@ *r) { - npy_intp n; - @ftype@ ar = npy_creal@c@(*a); - @ftype@ br = npy_creal@c@(*b); - @ftype@ ai = npy_cimag@c@(*a); - @ftype@ bi = npy_cimag@c@(*b); - - if (br == 0. && bi == 0.) { - *r = npy_cpack@c@(1., 0.); - return; - } - if (ar == 0. && ai == 0.) { - if (br > 0 && bi == 0) { - *r = npy_cpack@c@(0., 0.); - } - else { - volatile @ftype@ tmp = NPY_INFINITY; - /* NB: there are four complex zeros; c0 = (+-0, +-0), so that unlike - * for reals, c0**p, with `p` negative is in general - * ill-defined. - * - * c0**z with z complex is also ill-defined. - */ - *r = npy_cpack@c@(NPY_NAN, NPY_NAN); - - /* Raise invalid */ - tmp -= NPY_INFINITY; - ar = tmp; - } - return; - } - if (bi == 0 && (n=(npy_intp)br) == br) { - if (n == 1) { - /* unroll: handle inf better */ - *r = npy_cpack@c@(ar, ai); - return; - } - else if (n == 2) { - /* unroll: handle inf better */ - nc_prod@c@(a, a, r); - return; - } - else if (n == 3) { - /* unroll: handle inf better */ - nc_prod@c@(a, a, r); - nc_prod@c@(a, r, r); - return; - } - else if (n > -100 && n < 100) { - @ctype@ p, aa; - npy_intp mask = 1; - if (n < 0) n = -n; - aa = nc_1@c@; - p = npy_cpack@c@(ar, ai); - while (1) { - if (n & mask) - nc_prod@c@(&aa,&p,&aa); - mask <<= 1; - if (n < mask || mask <= 0) break; - nc_prod@c@(&p,&p,&p); - } - *r = npy_cpack@c@(npy_creal@c@(aa), npy_cimag@c@(aa)); - if (br < 0) nc_quot@c@(&nc_1@c@, r, r); - return; - } - } - - *r = npy_cpow@c@(*a, *b); + *r = npy_cpow@c@(*a, *b); return; } - -static void -nc_prodi@c@(@ctype@ *x, @ctype@ *r) -{ - @ftype@ xr = x->real; - r->real = -x->imag; - r->imag = xr; - return; -} - - static void nc_acos@c@(@ctype@ *x, @ctype@ *r) { - /* - * return nc_neg(nc_prodi(nc_log(nc_sum(x,nc_prod(nc_i, - * nc_sqrt(nc_diff(nc_1,nc_prod(x,x)))))))); - */ - nc_prod@c@(x,x,r); - nc_diff@c@(&nc_1@c@, r, r); - nc_sqrt@c@(r, r); - nc_prodi@c@(r, r); - nc_sum@c@(x, r, r); - nc_log@c@(r, r); - nc_prodi@c@(r, r); - nc_neg@c@(r, r); + *r = npy_cacos@c@(*x); return; } static void nc_acosh@c@(@ctype@ *x, @ctype@ *r) { - /* - * return nc_log(nc_sum(x, - * nc_prod(nc_sqrt(nc_sum(x,nc_1)), nc_sqrt(nc_diff(x,nc_1))))); - */ - @ctype@ t; - - nc_sum@c@(x, &nc_1@c@, &t); - nc_sqrt@c@(&t, &t); - nc_diff@c@(x, &nc_1@c@, r); - nc_sqrt@c@(r, r); - nc_prod@c@(&t, r, r); - nc_sum@c@(x, r, r); - nc_log@c@(r, r); + *r = npy_cacosh@c@(*x); return; } static void nc_asin@c@(@ctype@ *x, @ctype@ *r) { - /* - * return nc_neg(nc_prodi(nc_log(nc_sum(nc_prod(nc_i,x), - * nc_sqrt(nc_diff(nc_1,nc_prod(x,x))))))); - */ - if (fabs(x->real) > 1e-3 || fabs(x->imag) > 1e-3) { - @ctype@ a, *pa=&a; - nc_prod@c@(x, x, r); - nc_diff@c@(&nc_1@c@, r, r); - nc_sqrt@c@(r, r); - nc_prodi@c@(x, pa); - nc_sum@c@(pa, r, r); - nc_log@c@(r, r); - nc_prodi@c@(r, r); - nc_neg@c@(r, r); - } - else { - /* - * Small arguments: series expansion, to avoid loss of precision - * asin(x) = x [1 + (1/6) x^2 [1 + (9/20) x^2 [1 + ...]]] - * - * |x| < 1e-3 => |rel. error| < 1e-18 (f), 1e-24, 1e-36 (l) - */ - @ctype@ x2; - nc_prod@c@(x, x, &x2); - - *r = nc_1@c@; -#if @precision@ >= 3 - SERIES_HORNER_TERM@c@(r, &x2, 81.0@C@/110); - SERIES_HORNER_TERM@c@(r, &x2, 49.0@C@/72); -#endif -#if @precision@ >= 2 - SERIES_HORNER_TERM@c@(r, &x2, 25.0@C@/42); -#endif - SERIES_HORNER_TERM@c@(r, &x2, 9.0@C@/20); - SERIES_HORNER_TERM@c@(r, &x2, 1.0@C@/6); - nc_prod@c@(r, x, r); - } + *r = npy_casin@c@(*x); return; } @@ -472,134 +275,35 @@ nc_asin@c@(@ctype@ *x, @ctype@ *r) static void nc_asinh@c@(@ctype@ *x, @ctype@ *r) { - /* - * return nc_log(nc_sum(nc_sqrt(nc_sum(nc_1,nc_prod(x,x))),x)); - */ - if (fabs(x->real) > 1e-3 || fabs(x->imag) > 1e-3) { - nc_prod@c@(x, x, r); - nc_sum@c@(&nc_1@c@, r, r); - nc_sqrt@c@(r, r); - nc_sum@c@(r, x, r); - nc_log@c@(r, r); - } - else { - /* - * Small arguments: series expansion, to avoid loss of precision - * asinh(x) = x [1 - (1/6) x^2 [1 - (9/20) x^2 [1 - ...]]] - * - * |x| < 1e-3 => |rel. error| < 1e-18 (f), 1e-24, 1e-36 (l) - */ - @ctype@ x2; - nc_prod@c@(x, x, &x2); - - *r = nc_1@c@; -#if @precision@ >= 3 - SERIES_HORNER_TERM@c@(r, &x2, -81.0@C@/110); - SERIES_HORNER_TERM@c@(r, &x2, -49.0@C@/72); -#endif -#if @precision@ >= 2 - SERIES_HORNER_TERM@c@(r, &x2, -25.0@C@/42); -#endif - SERIES_HORNER_TERM@c@(r, &x2, -9.0@C@/20); - SERIES_HORNER_TERM@c@(r, &x2, -1.0@C@/6); - nc_prod@c@(r, x, r); - } + *r = npy_casinh@c@(*x); return; } static void nc_atan@c@(@ctype@ *x, @ctype@ *r) { - /* - * return nc_prod(nc_i2,nc_log(nc_quot(nc_sum(nc_i,x),nc_diff(nc_i,x)))); - */ - if (fabs(x->real) > 1e-3 || fabs(x->imag) > 1e-3) { - @ctype@ a, *pa=&a; - nc_diff@c@(&nc_i@c@, x, pa); - nc_sum@c@(&nc_i@c@, x, r); - nc_quot@c@(r, pa, r); - nc_log@c@(r,r); - nc_prod@c@(&nc_i2@c@, r, r); - } - else { - /* - * Small arguments: series expansion, to avoid loss of precision - * atan(x) = x [1 - (1/3) x^2 [1 - (3/5) x^2 [1 - ...]]] - * - * |x| < 1e-3 => |rel. error| < 1e-18 (f), 1e-24, 1e-36 (l) - */ - @ctype@ x2; - nc_prod@c@(x, x, &x2); - - *r = nc_1@c@; -#if @precision@ >= 3 - SERIES_HORNER_TERM@c@(r, &x2, -9.0@C@/11); - SERIES_HORNER_TERM@c@(r, &x2, -7.0@C@/9); -#endif -#if @precision@ >= 2 - SERIES_HORNER_TERM@c@(r, &x2, -5.0@C@/7); -#endif - SERIES_HORNER_TERM@c@(r, &x2, -3.0@C@/5); - SERIES_HORNER_TERM@c@(r, &x2, -1.0@C@/3); - nc_prod@c@(r, x, r); - } + *r = npy_catan@c@(*x); return; } static void nc_atanh@c@(@ctype@ *x, @ctype@ *r) { - /* - * return nc_prod(nc_half,nc_log(nc_quot(nc_sum(nc_1,x),nc_diff(nc_1,x)))); - */ - if (fabs(x->real) > 1e-3 || fabs(x->imag) > 1e-3) { - @ctype@ a, *pa=&a; - nc_diff@c@(&nc_1@c@, x, r); - nc_sum@c@(&nc_1@c@, x, pa); - nc_quot@c@(pa, r, r); - nc_log@c@(r, r); - nc_prod@c@(&nc_half@c@, r, r); - } - else { - /* - * Small arguments: series expansion, to avoid loss of precision - * atan(x) = x [1 + (1/3) x^2 [1 + (3/5) x^2 [1 + ...]]] - * - * |x| < 1e-3 => |rel. error| < 1e-18 (f), 1e-24, 1e-36 (l) - */ - @ctype@ x2; - nc_prod@c@(x, x, &x2); - - *r = nc_1@c@; -#if @precision@ >= 3 - SERIES_HORNER_TERM@c@(r, &x2, 9.0@C@/11); - SERIES_HORNER_TERM@c@(r, &x2, 7.0@C@/9); -#endif -#if @precision@ >= 2 - SERIES_HORNER_TERM@c@(r, &x2, 5.0@C@/7); -#endif - SERIES_HORNER_TERM@c@(r, &x2, 3.0@C@/5); - SERIES_HORNER_TERM@c@(r, &x2, 1.0@C@/3); - nc_prod@c@(r, x, r); - } + *r = npy_catanh@c@(*x); return; } static void nc_cos@c@(@ctype@ *x, @ctype@ *r) { - @ftype@ xr=x->real, xi=x->imag; - r->real = npy_cos@c@(xr)*npy_cosh@c@(xi); - r->imag = -npy_sin@c@(xr)*npy_sinh@c@(xi); + *r = npy_ccos@c@(*x); return; } static void nc_cosh@c@(@ctype@ *x, @ctype@ *r) { - @ftype@ xr=x->real, xi=x->imag; - r->real = npy_cos@c@(xi)*npy_cosh@c@(xr); - r->imag = npy_sin@c@(xi)*npy_sinh@c@(xr); + *r = npy_ccosh@c@(*x); return; } @@ -624,63 +328,29 @@ nc_log2@c@(@ctype@ *x, @ctype@ *r) static void nc_sin@c@(@ctype@ *x, @ctype@ *r) { - @ftype@ xr=x->real, xi=x->imag; - r->real = npy_sin@c@(xr)*npy_cosh@c@(xi); - r->imag = npy_cos@c@(xr)*npy_sinh@c@(xi); + *r = npy_csin@c@(*x); return; } static void nc_sinh@c@(@ctype@ *x, @ctype@ *r) { - @ftype@ xr=x->real, xi=x->imag; - r->real = npy_cos@c@(xi)*npy_sinh@c@(xr); - r->imag = npy_sin@c@(xi)*npy_cosh@c@(xr); + *r = npy_csinh@c@(*x); return; } static void nc_tan@c@(@ctype@ *x, @ctype@ *r) { - @ftype@ sr,cr,shi,chi; - @ftype@ rs,is,rc,ic; - @ftype@ d; - @ftype@ xr=x->real, xi=x->imag; - sr = npy_sin@c@(xr); - cr = npy_cos@c@(xr); - shi = npy_sinh@c@(xi); - chi = npy_cosh@c@(xi); - 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; - return; + *r = npy_ctan@c@(*x); + return; } static void nc_tanh@c@(@ctype@ *x, @ctype@ *r) { - @ftype@ si,ci,shr,chr; - @ftype@ rs,is,rc,ic; - @ftype@ d; - @ftype@ xr=x->real, xi=x->imag; - si = npy_sin@c@(xi); - ci = npy_cos@c@(xi); - shr = npy_sinh@c@(xr); - chr = npy_cosh@c@(xr); - 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; + *r = npy_ctanh@c@(*x); return; } -#undef SERIES_HORNER_TERM@c@ - /**end repeat**/ diff --git a/numpy/core/src/umath/ufunc_type_resolution.c b/numpy/core/src/umath/ufunc_type_resolution.c index ec28bb9e4..fe2e8cac3 100644 --- a/numpy/core/src/umath/ufunc_type_resolution.c +++ b/numpy/core/src/umath/ufunc_type_resolution.c @@ -1594,7 +1594,6 @@ set_ufunc_loop_data_types(PyUFuncObject *self, PyArrayObject **op, else if (op[i] != NULL && PyArray_DESCR(op[i])->type_num == type_nums[i]) { out_dtypes[i] = ensure_dtype_nbo(PyArray_DESCR(op[i])); - Py_XINCREF(out_dtypes[i]); } /* * For outputs, copy the dtype from op[0] if the type_num @@ -1603,7 +1602,6 @@ set_ufunc_loop_data_types(PyUFuncObject *self, PyArrayObject **op, else if (i >= nin && op[0] != NULL && PyArray_DESCR(op[0])->type_num == type_nums[i]) { out_dtypes[i] = ensure_dtype_nbo(PyArray_DESCR(op[0])); - Py_XINCREF(out_dtypes[i]); } /* Otherwise create a plain descr from the type number */ else { diff --git a/numpy/core/src/umath/umath_tests.c.src b/numpy/core/src/umath/umath_tests.c.src index 33d9846bd..509415711 100644 --- a/numpy/core/src/umath/umath_tests.c.src +++ b/numpy/core/src/umath/umath_tests.c.src @@ -202,12 +202,11 @@ static void INIT_OUTER_LOOP_2 npy_intp len_n = *dimensions++; npy_intp len_d = *dimensions++; - npy_intp len_p = *dimensions; npy_intp stride_n = *steps++; npy_intp stride_d = *steps++; npy_intp stride_p = *steps; - assert(len_n * (len_n - 1) / 2 == len_p); + assert(len_n * (len_n - 1) / 2 == *dimensions); BEGIN_OUTER_LOOP_2 const char *data_this = (const char *)args[0]; diff --git a/numpy/core/tests/test_numeric.py b/numpy/core/tests/test_numeric.py index b482bf230..9f0fb47e5 100644 --- a/numpy/core/tests/test_numeric.py +++ b/numpy/core/tests/test_numeric.py @@ -1574,6 +1574,11 @@ class TestAllclose(object): assert_(allclose(a, a)) + def test_equalnan(self): + x = np.array([1.0, np.nan]) + assert_(allclose(x, x, equal_nan=True)) + + class TestIsclose(object): rtol = 1e-5 atol = 1e-8 @@ -1664,17 +1669,25 @@ class TestIsclose(object): assert_array_equal(isclose(arr, arr, equal_nan=True), [True, True]) def test_masked_arrays(self): + # Make sure to test the output type when arguments are interchanged. + x = np.ma.masked_where([True, True, False], np.arange(3)) assert_(type(x) is type(isclose(2, x))) + assert_(type(x) is type(isclose(x, 2))) x = np.ma.masked_where([True, True, False], [nan, inf, nan]) assert_(type(x) is type(isclose(inf, x))) + assert_(type(x) is type(isclose(x, inf))) x = np.ma.masked_where([True, True, False], [nan, nan, nan]) y = isclose(nan, x, equal_nan=True) assert_(type(x) is type(y)) # Ensure that the mask isn't modified... assert_array_equal([True, True, False], y.mask) + y = isclose(x, nan, equal_nan=True) + assert_(type(x) is type(y)) + # Ensure that the mask isn't modified... + assert_array_equal([True, True, False], y.mask) x = np.ma.masked_where([True, True, False], [nan, nan, nan]) y = isclose(x, x, equal_nan=True) diff --git a/numpy/core/tests/test_records.py b/numpy/core/tests/test_records.py index 355e5480a..a7895a30a 100644 --- a/numpy/core/tests/test_records.py +++ b/numpy/core/tests/test_records.py @@ -73,10 +73,27 @@ class TestFromrecords(TestCase): assert_((mine.data2[i] == 0.0)) def test_recarray_from_repr(self): - x = np.rec.array([ (1, 2)], dtype=[('a', np.int8), ('b', np.int8)]) - y = eval("np." + repr(x)) - assert_(isinstance(y, np.recarray)) - assert_equal(y, x) + a = np.array([(1,'ABC'), (2, "DEF")], + dtype=[('foo', int), ('bar', 'S4')]) + recordarr = np.rec.array(a) + recarr = a.view(np.recarray) + recordview = a.view(np.dtype((np.record, a.dtype))) + + recordarr_r = eval("numpy." + repr(recordarr), {'numpy': np}) + recarr_r = eval("numpy." + repr(recarr), {'numpy': np}) + recordview_r = eval("numpy." + repr(recordview), {'numpy': np}) + + assert_equal(type(recordarr_r), np.recarray) + assert_equal(recordarr_r.dtype.type, np.record) + assert_equal(recordarr, recordarr_r) + + assert_equal(type(recarr_r), np.recarray) + assert_equal(recarr_r.dtype.type, np.void) + assert_equal(recarr, recarr_r) + + assert_equal(type(recordview_r), np.ndarray) + assert_equal(recordview.dtype.type, np.record) + assert_equal(recordview, recordview_r) def test_recarray_from_names(self): ra = np.rec.array([ @@ -124,6 +141,28 @@ class TestFromrecords(TestCase): assert_equal(a.b, ['a', 'bbb']) assert_equal(a[-1].b, 'bbb') + def test_recarray_stringtypes(self): + # Issue #3993 + a = np.array([('abc ', 1), ('abc', 2)], + dtype=[('foo', 'S4'), ('bar', int)]) + a = a.view(np.recarray) + assert_equal(a.foo[0] == a.foo[1], False) + + def test_recarray_returntypes(self): + a = np.rec.array([('abc ', (1,1), 1), ('abc', (2,3), 1)], + dtype=[('foo', 'S4'), + ('bar', [('A', int), ('B', int)]), + ('baz', int)]) + assert_equal(type(a.foo), np.ndarray) + assert_equal(type(a['foo']), np.ndarray) + assert_equal(type(a.bar), np.recarray) + assert_equal(type(a['bar']), np.recarray) + assert_equal(a.bar.dtype.type, np.record) + assert_equal(type(a.baz), np.ndarray) + assert_equal(type(a['baz']), np.ndarray) + assert_equal(type(a[0].bar), np.record) + assert_equal(a[0].bar.A, 1) + class TestRecord(TestCase): def setUp(self): diff --git a/numpy/core/tests/test_umath.py b/numpy/core/tests/test_umath.py index c71b7b658..092953872 100644 --- a/numpy/core/tests/test_umath.py +++ b/numpy/core/tests/test_umath.py @@ -1363,45 +1363,53 @@ class TestComplexFunctions(object): def test_branch_cuts(self): # check branch cuts and continuity on them - yield _check_branch_cut, np.log, -0.5, 1j, 1, -1 - yield _check_branch_cut, np.log2, -0.5, 1j, 1, -1 - yield _check_branch_cut, np.log10, -0.5, 1j, 1, -1 - yield _check_branch_cut, np.log1p, -1.5, 1j, 1, -1 - yield _check_branch_cut, np.sqrt, -0.5, 1j, 1, -1 + yield _check_branch_cut, np.log, -0.5, 1j, 1, -1, True + yield _check_branch_cut, np.log2, -0.5, 1j, 1, -1, True + yield _check_branch_cut, np.log10, -0.5, 1j, 1, -1, True + yield _check_branch_cut, np.log1p, -1.5, 1j, 1, -1, True + yield _check_branch_cut, np.sqrt, -0.5, 1j, 1, -1, True - yield _check_branch_cut, np.arcsin, [ -2, 2], [1j, -1j], 1, -1 - yield _check_branch_cut, np.arccos, [ -2, 2], [1j, -1j], 1, -1 - yield _check_branch_cut, np.arctan, [-2j, 2j], [1, -1 ], -1, 1 + yield _check_branch_cut, np.arcsin, [ -2, 2], [1j, 1j], 1, -1, True + yield _check_branch_cut, np.arccos, [ -2, 2], [1j, 1j], 1, -1, True + yield _check_branch_cut, np.arctan, [0-2j, 2j], [1, 1 ], -1, 1, True - yield _check_branch_cut, np.arcsinh, [-2j, 2j], [-1, 1], -1, 1 - yield _check_branch_cut, np.arccosh, [ -1, 0.5], [1j, 1j], 1, -1 - yield _check_branch_cut, np.arctanh, [ -2, 2], [1j, -1j], 1, -1 + yield _check_branch_cut, np.arcsinh, [0-2j, 2j], [1, 1], -1, 1, True + yield _check_branch_cut, np.arccosh, [ -1, 0.5], [1j, 1j], 1, -1, True + yield _check_branch_cut, np.arctanh, [ -2, 2], [1j, 1j], 1, -1, True # check against bogus branch cuts: assert continuity between quadrants - yield _check_branch_cut, np.arcsin, [-2j, 2j], [ 1, 1], 1, 1 - yield _check_branch_cut, np.arccos, [-2j, 2j], [ 1, 1], 1, 1 + yield _check_branch_cut, np.arcsin, [0-2j, 2j], [ 1, 1], 1, 1 + yield _check_branch_cut, np.arccos, [0-2j, 2j], [ 1, 1], 1, 1 yield _check_branch_cut, np.arctan, [ -2, 2], [1j, 1j], 1, 1 yield _check_branch_cut, np.arcsinh, [ -2, 2, 0], [1j, 1j, 1 ], 1, 1 - yield _check_branch_cut, np.arccosh, [-2j, 2j, 2], [1, 1, 1j], 1, 1 - yield _check_branch_cut, np.arctanh, [-2j, 2j, 0], [1, 1, 1j], 1, 1 + yield _check_branch_cut, np.arccosh, [0-2j, 2j, 2], [1, 1, 1j], 1, 1 + yield _check_branch_cut, np.arctanh, [0-2j, 2j, 0], [1, 1, 1j], 1, 1 - @dec.knownfailureif(True, "These branch cuts are known to fail") - def test_branch_cuts_failing(self): - # XXX: signed zero not OK with ICC on 64-bit platform for log, see - # http://permalink.gmane.org/gmane.comp.python.numeric.general/25335 - yield _check_branch_cut, np.log, -0.5, 1j, 1, -1, True - yield _check_branch_cut, np.log2, -0.5, 1j, 1, -1, True - yield _check_branch_cut, np.log10, -0.5, 1j, 1, -1, True - yield _check_branch_cut, np.log1p, -1.5, 1j, 1, -1, True - # XXX: signed zeros are not OK for sqrt or for the arc* functions - yield _check_branch_cut, np.sqrt, -0.5, 1j, 1, -1, True - yield _check_branch_cut, np.arcsin, [ -2, 2], [1j, -1j], 1, -1, True - yield _check_branch_cut, np.arccos, [ -2, 2], [1j, -1j], 1, -1, True - yield _check_branch_cut, np.arctan, [-2j, 2j], [1, -1 ], -1, 1, True - yield _check_branch_cut, np.arcsinh, [-2j, 2j], [-1, 1], -1, 1, True - yield _check_branch_cut, np.arccosh, [ -1, 0.5], [1j, 1j], 1, -1, True - yield _check_branch_cut, np.arctanh, [ -2, 2], [1j, -1j], 1, -1, True + def test_branch_cuts_complex64(self): + # check branch cuts and continuity on them + yield _check_branch_cut, np.log, -0.5, 1j, 1, -1, True, np.complex64 + yield _check_branch_cut, np.log2, -0.5, 1j, 1, -1, True, np.complex64 + yield _check_branch_cut, np.log10, -0.5, 1j, 1, -1, True, np.complex64 + yield _check_branch_cut, np.log1p, -1.5, 1j, 1, -1, True, np.complex64 + yield _check_branch_cut, np.sqrt, -0.5, 1j, 1, -1, True, np.complex64 + + yield _check_branch_cut, np.arcsin, [ -2, 2], [1j, 1j], 1, -1, True, np.complex64 + yield _check_branch_cut, np.arccos, [ -2, 2], [1j, 1j], 1, -1, True, np.complex64 + yield _check_branch_cut, np.arctan, [0-2j, 2j], [1, 1 ], -1, 1, True, np.complex64 + + yield _check_branch_cut, np.arcsinh, [0-2j, 2j], [1, 1], -1, 1, True, np.complex64 + yield _check_branch_cut, np.arccosh, [ -1, 0.5], [1j, 1j], 1, -1, True, np.complex64 + yield _check_branch_cut, np.arctanh, [ -2, 2], [1j, 1j], 1, -1, True, np.complex64 + + # check against bogus branch cuts: assert continuity between quadrants + yield _check_branch_cut, np.arcsin, [0-2j, 2j], [ 1, 1], 1, 1, False, np.complex64 + yield _check_branch_cut, np.arccos, [0-2j, 2j], [ 1, 1], 1, 1, False, np.complex64 + yield _check_branch_cut, np.arctan, [ -2, 2], [1j, 1j], 1, 1, False, np.complex64 + + yield _check_branch_cut, np.arcsinh, [ -2, 2, 0], [1j, 1j, 1 ], 1, 1, False, np.complex64 + yield _check_branch_cut, np.arccosh, [0-2j, 2j, 2], [1, 1, 1j], 1, 1, False, np.complex64 + yield _check_branch_cut, np.arctanh, [0-2j, 2j, 0], [1, 1, 1j], 1, 1, False, np.complex64 def test_against_cmath(self): import cmath, sys @@ -1466,7 +1474,7 @@ class TestComplexFunctions(object): # So, give more leeway for long complex tests here: check(x_series, 50*eps) else: - check(x_series, 2*eps) + check(x_series, 2.1*eps) check(x_basic, 2*eps/1e-3) # Check a few points @@ -1565,8 +1573,12 @@ def _check_branch_cut(f, x0, dx, re_sign=1, im_sign=-1, sig_zero_ok=False, x0 = np.atleast_1d(x0).astype(dtype) dx = np.atleast_1d(dx).astype(dtype) - scale = np.finfo(dtype).eps * 1e3 - atol = 1e-4 + if np.dtype(dtype).char == 'F': + scale = np.finfo(dtype).eps * 1e2 + atol = np.float32(1e-2) + else: + scale = np.finfo(dtype).eps * 1e3 + atol = 1e-4 y0 = f(x0) yp = f(x0 + dx*scale*np.absolute(x0)/np.absolute(dx)) @@ -1581,16 +1593,20 @@ def _check_branch_cut(f, x0, dx, re_sign=1, im_sign=-1, sig_zero_ok=False, # check that signed zeros also work as a displacement jr = (x0.real == 0) & (dx.real != 0) ji = (x0.imag == 0) & (dx.imag != 0) - - x = -x0 - x.real[jr] = 0.*dx.real - x.imag[ji] = 0.*dx.imag - x = -x - ym = f(x) - ym = ym[jr | ji] - y0 = y0[jr | ji] - assert_(np.all(np.absolute(y0.real - ym.real*re_sign) < atol), (y0, ym)) - assert_(np.all(np.absolute(y0.imag - ym.imag*im_sign) < atol), (y0, ym)) + if np.any(jr): + x = x0[jr] + x.real = np.NZERO + ym = f(x) + assert_(np.all(np.absolute(y0[jr].real - ym.real*re_sign) < atol), (y0[jr], ym)) + assert_(np.all(np.absolute(y0[jr].imag - ym.imag*im_sign) < atol), (y0[jr], ym)) + + + if np.any(ji): + x = x0[ji] + x.imag = np.NZERO + ym = f(x) + assert_(np.all(np.absolute(y0[ji].real - ym.real*re_sign) < atol), (y0[ji], ym)) + assert_(np.all(np.absolute(y0[ji].imag - ym.imag*im_sign) < atol), (y0[ji], ym)) def test_copysign(): assert_(np.copysign(1, -1) == -1) diff --git a/numpy/distutils/fcompiler/intel.py b/numpy/distutils/fcompiler/intel.py index a80e525e3..f76174c7a 100644 --- a/numpy/distutils/fcompiler/intel.py +++ b/numpy/distutils/fcompiler/intel.py @@ -152,7 +152,7 @@ class IntelVisualFCompiler(BaseIntelFCompiler): module_include_switch = '/I' def get_flags(self): - opt = ['/nologo', '/MD', '/nbs', '/Qlowercase', '/us'] + opt = ['/nologo', '/MD', '/nbs', '/names:lowercase', '/assume:underscore'] return opt def get_flags_free(self): diff --git a/numpy/distutils/lib2def.py b/numpy/distutils/lib2def.py index 7316547a3..0a5364566 100644 --- a/numpy/distutils/lib2def.py +++ b/numpy/distutils/lib2def.py @@ -66,7 +66,7 @@ def getnm(nm_cmd = ['nm', '-Cs', 'python%s.lib' % py_ver]): """Returns the output of nm_cmd via a pipe. nm_output = getnam(nm_cmd = 'nm -Cs py_lib')""" - f = subprocess.Popen(nm_cmd, shell=True, stdout=subprocess.PIPE) + f = subprocess.Popen(nm_cmd, shell=True, stdout=subprocess.PIPE, universal_newlines=True) nm_output = f.stdout.read() f.stdout.close() return nm_output diff --git a/numpy/doc/structured_arrays.py b/numpy/doc/structured_arrays.py index f2329827e..73bf3b317 100644 --- a/numpy/doc/structured_arrays.py +++ b/numpy/doc/structured_arrays.py @@ -268,6 +268,10 @@ array if the field has a structured type but as a plain ndarray otherwise. :: >>> type(recordarr.bar) <class 'numpy.core.records.recarray'> +Note that if a field has the same name as an ndarray attribute, the ndarray +attribute takes precedence. Such fields will be inaccessible by attribute but +may still be accessed by index. + Partial Attribute Access ------------------------ @@ -288,7 +292,7 @@ such a view do not have field attributes:: To use the np.record dtype only, convert the dtype using the (base_class, dtype) form described in numpy.dtype. This type of view is rarely used. :: - >>> arr_records = arr.view(dtype(np.record, arr.dtype)) + >>> arr_records = arr.view(dtype((np.record, arr.dtype))) In documentation, the term 'structured array' will refer to objects of type np.ndarray with structured dtype, 'record array' will refer to structured diff --git a/numpy/fft/fftpack.c b/numpy/fft/fftpack.c index 355076a30..277f49f07 100644 --- a/numpy/fft/fftpack.c +++ b/numpy/fft/fftpack.c @@ -1,18 +1,15 @@ /* -fftpack.c : A set of FFT routines in C. -Algorithmically based on Fortran-77 FFTPACK by Paul N. Swarztrauber (Version 4, 1985). - + * fftpack.c : A set of FFT routines in C. + * Algorithmically based on Fortran-77 FFTPACK by Paul N. Swarztrauber (Version 4, 1985). */ +#define NPY_NO_DEPRECATED_API NPY_API_VERSION -/* isign is +1 for backward and -1 for forward transforms */ - +#include <Python.h> #include <math.h> #include <stdio.h> -#include <Python.h> #include <numpy/ndarraytypes.h> #define DOUBLE - #ifdef DOUBLE #define Treal double #else diff --git a/numpy/fft/fftpack_litemodule.c b/numpy/fft/fftpack_litemodule.c index 7b37ce9b7..e895d0efe 100644 --- a/numpy/fft/fftpack_litemodule.c +++ b/numpy/fft/fftpack_litemodule.c @@ -149,7 +149,7 @@ fftpack_rfftf(PyObject *NPY_UNUSED(self), PyObject *args) PyObject *op1, *op2; PyArrayObject *data, *ret; PyArray_Descr *descr; - double *wsave, *dptr, *rptr; + double *wsave = NULL, *dptr, *rptr; npy_intp nsave; int npts, nrepeats, i, rstep; @@ -166,6 +166,9 @@ fftpack_rfftf(PyObject *NPY_UNUSED(self), PyObject *args) PyArray_DIMS(data)[PyArray_NDIM(data) - 1] = npts/2 + 1; ret = (PyArrayObject *)PyArray_Zeros(PyArray_NDIM(data), PyArray_DIMS(data), PyArray_DescrFromType(NPY_CDOUBLE), 0); + if (ret == NULL) { + goto fail; + } PyArray_DIMS(data)[PyArray_NDIM(data) - 1] = npts; rstep = PyArray_DIM(ret, PyArray_NDIM(ret) - 1)*2; diff --git a/numpy/lib/_iotools.py b/numpy/lib/_iotools.py index f2adcda10..316704b42 100644 --- a/numpy/lib/_iotools.py +++ b/numpy/lib/_iotools.py @@ -320,12 +320,13 @@ class NameValidator(object): # Process the case option ..... if (case_sensitive is None) or (case_sensitive is True): self.case_converter = lambda x: x - elif (case_sensitive is False) or ('u' in case_sensitive): + elif (case_sensitive is False) or case_sensitive.startswith('u'): self.case_converter = lambda x: x.upper() - elif 'l' in case_sensitive: + elif case_sensitive.startswith('l'): self.case_converter = lambda x: x.lower() else: - self.case_converter = lambda x: x + msg = 'unrecognized case_sensitive value %s.' % case_sensitive + raise ValueError(msg) # self.replace_space = replace_space diff --git a/numpy/lib/arraysetops.py b/numpy/lib/arraysetops.py index d3b6119f4..cb24eb24e 100644 --- a/numpy/lib/arraysetops.py +++ b/numpy/lib/arraysetops.py @@ -392,12 +392,13 @@ def in1d(ar1, ar2, assume_unique=False, invert=False): else: bool_ar = (sar[1:] == sar[:-1]) flag = np.concatenate((bool_ar, [invert])) - indx = order.argsort(kind='mergesort')[:len(ar1)] + ret = np.empty(ar.shape, dtype=bool) + ret[order] = flag if assume_unique: - return flag[indx] + return ret[:len(ar1)] else: - return flag[indx][rev_idx] + return ret[rev_idx] def union1d(ar1, ar2): """ diff --git a/numpy/lib/financial.py b/numpy/lib/financial.py index baff8b0b6..a7e4e60b6 100644 --- a/numpy/lib/financial.py +++ b/numpy/lib/financial.py @@ -208,11 +208,11 @@ def pmt(rate, nper, pv, fv=0, when='end'): """ when = _convert_when(when) (rate, nper, pv, fv, when) = map(np.asarray, [rate, nper, pv, fv, when]) - temp = (1+rate)**nper - miter = np.broadcast(rate, nper, pv, fv, when) - zer = np.zeros(miter.shape) - fact = np.where(rate == zer, nper + zer, - (1 + rate*when)*(temp - 1)/rate + zer) + temp = (1 + rate)**nper + mask = (rate == 0.0) + np.copyto(rate, 1.0, where=mask) + z = np.zeros(np.broadcast(rate, nper, pv, fv, when).shape) + fact = np.where(mask != z, nper + z, (1 + rate*when)*(temp - 1)/rate + z) return -(fv + pv*temp) / fact def nper(rate, pmt, pv, fv=0, when='end'): diff --git a/numpy/lib/npyio.py b/numpy/lib/npyio.py index 9d539d6ac..0632ba1f8 100644 --- a/numpy/lib/npyio.py +++ b/numpy/lib/npyio.py @@ -621,6 +621,13 @@ def _savez(file, args, kwds, compress): def _getconv(dtype): """ Find the correct dtype converter. Adapted from matplotlib """ + + def floatconv(x): + x.lower() + if b'0x' in x: + return float.fromhex(asstr(x)) + return float(x) + typ = dtype.type if issubclass(typ, np.bool_): return lambda x: bool(int(x)) @@ -631,7 +638,7 @@ def _getconv(dtype): if issubclass(typ, np.integer): return lambda x: int(float(x)) elif issubclass(typ, np.floating): - return float + return floatconv elif issubclass(typ, np.complex): return complex elif issubclass(typ, np.bytes_): @@ -706,6 +713,10 @@ def loadtxt(fname, dtype=float, comments='#', delimiter=None, `genfromtxt` function provides more sophisticated handling of, e.g., lines with missing values. + .. versionadded:: 1.10.0 + The strings produced by the Python float.hex method can be used as + input for floats. + Examples -------- >>> from StringIO import StringIO # StringIO behaves like a file object @@ -1204,7 +1215,8 @@ def genfromtxt(fname, dtype=float, comments='#', delimiter=None, usecols=None, names=None, excludelist=None, deletechars=None, replace_space='_', autostrip=False, case_sensitive=True, defaultfmt="f%i", - unpack=None, usemask=False, loose=True, invalid_raise=True): + unpack=None, usemask=False, loose=True, invalid_raise=True, + max_rows=None): """ Load data from a text file, with missing values handled as specified. @@ -1285,6 +1297,12 @@ def genfromtxt(fname, dtype=float, comments='#', delimiter=None, If True, an exception is raised if an inconsistency is detected in the number of columns. If False, a warning is emitted and the offending lines are skipped. + max_rows : int, optional + The maximum number of rows to read. Must not be used with skip_footer + at the same time. If given, the value must be at least 1. Default is + to read the entire file. + + .. versionadded:: 1.10.0 Returns ------- @@ -1353,6 +1371,14 @@ def genfromtxt(fname, dtype=float, comments='#', delimiter=None, dtype=[('intvar', '<i8'), ('fltvar', '<f8'), ('strvar', '|S5')]) """ + if max_rows is not None: + if skip_footer: + raise ValueError( + "The keywords 'skip_footer' and 'max_rows' can not be " + "specified at the same time.") + if max_rows < 1: + raise ValueError("'max_rows' must be at least 1.") + # Py3 data conversions to bytes, for convenience if comments is not None: comments = asbytes(comments) @@ -1451,7 +1477,11 @@ def genfromtxt(fname, dtype=float, comments='#', delimiter=None, names = validate_names(names) # Get the dtype if dtype is not None: - dtype = easy_dtype(dtype, defaultfmt=defaultfmt, names=names) + dtype = easy_dtype(dtype, defaultfmt=defaultfmt, names=names, + excludelist=excludelist, + deletechars=deletechars, + case_sensitive=case_sensitive, + replace_space=replace_space) # Make sure the names is a list (for 2.5) if names is not None: names = list(names) @@ -1647,8 +1677,8 @@ def genfromtxt(fname, dtype=float, comments='#', delimiter=None, # Skip an empty line if nbvalues == 0: continue - # Select only the columns we need if usecols: + # Select only the columns we need try: values = [values[_] for _ in usecols] except IndexError: @@ -1661,7 +1691,10 @@ def genfromtxt(fname, dtype=float, comments='#', delimiter=None, append_to_rows(tuple(values)) if usemask: append_to_masks(tuple([v.strip() in m - for (v, m) in zip(values, missing_values)])) + for (v, m) in zip(values, + missing_values)])) + if len(rows) == max_rows: + break if own_fhd: fhd.close() diff --git a/numpy/lib/setup.py b/numpy/lib/setup.py new file mode 100644 index 000000000..d342410b8 --- /dev/null +++ b/numpy/lib/setup.py @@ -0,0 +1,12 @@ +from __future__ import division, print_function + +def configuration(parent_package='',top_path=None): + from numpy.distutils.misc_util import Configuration + + config = Configuration('lib', parent_package, top_path) + config.add_data_dir('tests') + return config + +if __name__ == '__main__': + from numpy.distutils.core import setup + setup(configuration=configuration) diff --git a/numpy/lib/shape_base.py b/numpy/lib/shape_base.py index 70fa3ab03..2d18c5bc8 100644 --- a/numpy/lib/shape_base.py +++ b/numpy/lib/shape_base.py @@ -711,7 +711,7 @@ def kron(a, b): Notes ----- - The function assumes that the number of dimenensions of `a` and `b` + The function assumes that the number of dimensions of `a` and `b` are the same, if necessary prepending the smallest with ones. If `a.shape = (r0,r1,..,rN)` and `b.shape = (s0,s1,...,sN)`, the Kronecker product has shape `(r0*s0, r1*s1, ..., rN*SN)`. diff --git a/numpy/lib/tests/test__iotools.py b/numpy/lib/tests/test__iotools.py index 92ca1c973..060f815d5 100644 --- a/numpy/lib/tests/test__iotools.py +++ b/numpy/lib/tests/test__iotools.py @@ -7,7 +7,8 @@ from datetime import date import numpy as np from numpy.compat import asbytes, asbytes_nested from numpy.testing import ( - run_module_suite, TestCase, assert_, assert_equal, assert_allclose + run_module_suite, TestCase, assert_, assert_equal, assert_allclose, + assert_raises ) from numpy.lib._iotools import ( LineSplitter, NameValidator, StringConverter, @@ -93,6 +94,9 @@ class TestNameValidator(TestCase): test = NameValidator(case_sensitive='lower').validate(names) assert_equal(test, ['a', 'a_1', 'b', 'c']) + # check exceptions + assert_raises(ValueError, NameValidator, case_sensitive='foobar') + def test_excludelist(self): "Test excludelist" names = ['dates', 'data', 'Other Data', 'mask'] diff --git a/numpy/lib/tests/test_financial.py b/numpy/lib/tests/test_financial.py index a4b9cfe2e..baa785424 100644 --- a/numpy/lib/tests/test_financial.py +++ b/numpy/lib/tests/test_financial.py @@ -2,7 +2,8 @@ from __future__ import division, absolute_import, print_function import numpy as np from numpy.testing import ( - run_module_suite, TestCase, assert_, assert_almost_equal + run_module_suite, TestCase, assert_, assert_almost_equal, + assert_allclose ) @@ -13,35 +14,37 @@ class TestFinancial(TestCase): def test_irr(self): v = [-150000, 15000, 25000, 35000, 45000, 60000] - assert_almost_equal(np.irr(v), - 0.0524, 2) + assert_almost_equal(np.irr(v), 0.0524, 2) v = [-100, 0, 0, 74] - assert_almost_equal(np.irr(v), - -0.0955, 2) + assert_almost_equal(np.irr(v), -0.0955, 2) v = [-100, 39, 59, 55, 20] - assert_almost_equal(np.irr(v), - 0.28095, 2) + assert_almost_equal(np.irr(v), 0.28095, 2) v = [-100, 100, 0, -7] - assert_almost_equal(np.irr(v), - -0.0833, 2) + assert_almost_equal(np.irr(v), -0.0833, 2) v = [-100, 100, 0, 7] - assert_almost_equal(np.irr(v), - 0.06206, 2) + assert_almost_equal(np.irr(v), 0.06206, 2) v = [-5, 10.5, 1, -8, 1] - assert_almost_equal(np.irr(v), - 0.0886, 2) + assert_almost_equal(np.irr(v), 0.0886, 2) def test_pv(self): - assert_almost_equal(np.pv(0.07, 20, 12000, 0), - -127128.17, 2) + assert_almost_equal(np.pv(0.07, 20, 12000, 0), -127128.17, 2) def test_fv(self): - assert_almost_equal(np.fv(0.075, 20, -2000, 0, 0), - 86609.36, 2) + assert_almost_equal(np.fv(0.075, 20, -2000, 0, 0), 86609.36, 2) def test_pmt(self): - assert_almost_equal(np.pmt(0.08/12, 5*12, 15000), - -304.146, 3) + res = np.pmt(0.08/12, 5*12, 15000) + tgt = -304.145914 + assert_allclose(res, tgt) + # Test the edge case where rate == 0.0 + res = np.pmt(0.0, 5*12, 15000) + tgt = -250.0 + assert_allclose(res, tgt) + # Test the case where we use broadcast and + # the arguments passed in are arrays. + res = np.pmt([[0.0, 0.8],[0.3, 0.8]],[12, 3],[2000, 20000]) + tgt = np.array([[-166.66667, -19311.258],[-626.90814, -19311.258]]) + assert_allclose(res, tgt) def test_ppmt(self): np.round(np.ppmt(0.1/12, 1, 60, 55000), 2) == 710.25 diff --git a/numpy/lib/tests/test_io.py b/numpy/lib/tests/test_io.py index 81bddfadd..7054ab1fe 100644 --- a/numpy/lib/tests/test_io.py +++ b/numpy/lib/tests/test_io.py @@ -694,6 +694,19 @@ class TestLoadTxt(TestCase): res = np.loadtxt(c, dtype=np.int64) assert_equal(res, tgt) + def test_from_float_hex(self): + # IEEE doubles and floats only, otherwise the float32 + # conversion may fail. + tgt = np.logspace(-10, 10, 5).astype(np.float32) + tgt = np.hstack((tgt, -tgt)).astype(np.float) + inp = '\n'.join(map(float.hex, tgt)) + c = TextIO() + c.write(inp) + for dt in [np.float, np.float32]: + c.seek(0) + res = np.loadtxt(c, dtype=dt) + assert_equal(res, tgt, err_msg="%s" % dt) + def test_universal_newline(self): f, name = mkstemp() os.write(f, b'1 21\r3 42\r') @@ -1107,13 +1120,13 @@ M 33 21.99 def test_dtype_with_converters_and_usecols(self): dstr = "1,5,-1,1:1\n2,8,-1,1:n\n3,3,-2,m:n\n" dmap = {'1:1':0, '1:n':1, 'm:1':2, 'm:n':3} - dtyp = [('E1','i4'),('E2','i4'),('E3','i2'),('N', 'i1')] + dtyp = [('e1','i4'),('e2','i4'),('e3','i2'),('n', 'i1')] conv = {0: int, 1: int, 2: int, 3: lambda r: dmap[r.decode()]} test = np.recfromcsv(TextIO(dstr,), dtype=dtyp, delimiter=',', names=None, converters=conv) control = np.rec.array([[1,5,-1,0], [2,8,-1,1], [3,3,-2,3]], dtype=dtyp) assert_equal(test, control) - dtyp = [('E1','i4'),('E2','i4'),('N', 'i1')] + dtyp = [('e1','i4'),('e2','i4'),('n', 'i1')] test = np.recfromcsv(TextIO(dstr,), dtype=dtyp, delimiter=',', usecols=(0,1,3), names=None, converters=conv) control = np.rec.array([[1,5,0], [2,8,1], [3,3,3]], dtype=dtyp) @@ -1514,6 +1527,30 @@ M 33 21.99 ctrl = np.array((1, 2, 3.14), dtype=ctrl_dtype) assert_equal(test, ctrl) + def test_replace_space_known_dtype(self): + "Test the 'replace_space' (and related) options when dtype != None" + txt = "A.A, B (B), C:C\n1, 2, 3" + # Test default: replace ' ' by '_' and delete non-alphanum chars + test = np.genfromtxt(TextIO(txt), + delimiter=",", names=True, dtype=int) + ctrl_dtype = [("AA", int), ("B_B", int), ("CC", int)] + ctrl = np.array((1, 2, 3), dtype=ctrl_dtype) + assert_equal(test, ctrl) + # Test: no replace, no delete + test = np.genfromtxt(TextIO(txt), + delimiter=",", names=True, dtype=int, + replace_space='', deletechars='') + ctrl_dtype = [("A.A", int), ("B (B)", int), ("C:C", int)] + ctrl = np.array((1, 2, 3), dtype=ctrl_dtype) + assert_equal(test, ctrl) + # Test: no delete (spaces are replaced by _) + test = np.genfromtxt(TextIO(txt), + delimiter=",", names=True, dtype=int, + deletechars='') + ctrl_dtype = [("A.A", int), ("B_(B)", int), ("C:C", int)] + ctrl = np.array((1, 2, 3), dtype=ctrl_dtype) + assert_equal(test, ctrl) + def test_incomplete_names(self): "Test w/ incomplete names" data = "A,,C\n0,1,2\n3,4,5" @@ -1641,6 +1678,60 @@ M 33 21.99 self.assertTrue(isinstance(test, np.recarray)) assert_equal(test, control) + def test_max_rows(self): + # Test the `max_rows` keyword argument. + data = '1 2\n3 4\n5 6\n7 8\n9 10\n' + txt = TextIO(data) + a1 = np.genfromtxt(txt, max_rows=3) + a2 = np.genfromtxt(txt) + assert_equal(a1, [[1, 2], [3, 4], [5, 6]]) + assert_equal(a2, [[7, 8], [9, 10]]) + + # max_rows must be at least 1. + assert_raises(ValueError, np.genfromtxt, TextIO(data), max_rows=0) + + # An input with several invalid rows. + data = '1 1\n2 2\n0 \n3 3\n4 4\n5 \n6 \n7 \n' + + test = np.genfromtxt(TextIO(data), max_rows=2) + control = np.array([[1., 1.], [2., 2.]]) + assert_equal(test, control) + + # Test keywords conflict + assert_raises(ValueError, np.genfromtxt, TextIO(data), skip_footer=1, + max_rows=4) + + # Test with invalid value + assert_raises(ValueError, np.genfromtxt, TextIO(data), max_rows=4) + + # Test with invalid not raise + with warnings.catch_warnings(): + warnings.filterwarnings("ignore") + + test = np.genfromtxt(TextIO(data), max_rows=4, invalid_raise=False) + control = np.array([[1., 1.], [2., 2.], [3., 3.], [4., 4.]]) + assert_equal(test, control) + + test = np.genfromtxt(TextIO(data), max_rows=5, invalid_raise=False) + control = np.array([[1., 1.], [2., 2.], [3., 3.], [4., 4.]]) + assert_equal(test, control) + + # Structured array with field names. + data = 'a b\n#c d\n1 1\n2 2\n#0 \n3 3\n4 4\n5 5\n' + + # Test with header, names and comments + txt = TextIO(data) + test = np.genfromtxt(txt, skip_header=1, max_rows=3, names=True) + control = np.array([(1.0, 1.0), (2.0, 2.0), (3.0, 3.0)], + dtype=[('c', '<f8'), ('d', '<f8')]) + assert_equal(test, control) + # To continue reading the same "file", don't use skip_header or + # names, and use the previously determined dtype. + test = np.genfromtxt(txt, max_rows=None, dtype=test.dtype) + control = np.array([(4.0, 4.0), (5.0, 5.0)], + dtype=[('c', '<f8'), ('d', '<f8')]) + assert_equal(test, control) + def test_gft_using_filename(self): # Test that we can load data from a filename as well as a file object wanted = np.arange(6).reshape((2, 3)) diff --git a/numpy/lib/tests/test_twodim_base.py b/numpy/lib/tests/test_twodim_base.py index 739061a5d..b65a8df97 100644 --- a/numpy/lib/tests/test_twodim_base.py +++ b/numpy/lib/tests/test_twodim_base.py @@ -265,6 +265,37 @@ class TestHistogram2d(TestCase): a, edge1, edge2 = histogram2d([], [], bins=4) assert_array_max_ulp(a, np.zeros((4, 4))) + def test_binparameter_combination(self): + x = array( + [0, 0.09207008, 0.64575234, 0.12875982, 0.47390599, + 0.59944483, 1]) + y = array( + [0, 0.14344267, 0.48988575, 0.30558665, 0.44700682, + 0.15886423, 1]) + edges = (0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1) + H, xe, ye = histogram2d(x, y, (edges, 4)) + answer = array( + [[ 2., 0., 0., 0.], + [ 0., 1., 0., 0.], + [ 0., 0., 0., 0.], + [ 0., 0., 0., 0.], + [ 0., 1., 0., 0.], + [ 1., 0., 0., 0.], + [ 0., 1., 0., 0.], + [ 0., 0., 0., 0.], + [ 0., 0., 0., 0.], + [ 0., 0., 0., 1.]]) + assert_array_equal(H, answer) + assert_array_equal(ye, array([0., 0.25, 0.5, 0.75, 1])) + H, xe, ye = histogram2d(x, y, (4, edges)) + answer = array( + [[ 1., 1., 0., 1., 0., 0., 0., 0., 0., 0.], + [ 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.], + [ 0., 1., 0., 0., 1., 0., 0., 0., 0., 0.], + [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.]]) + assert_array_equal(H, answer) + assert_array_equal(xe, array([0., 0.25, 0.5, 0.75, 1])) + class TestTri(TestCase): def test_dtype(self): diff --git a/numpy/lib/twodim_base.py b/numpy/lib/twodim_base.py index 24dc3cced..464ffd914 100644 --- a/numpy/lib/twodim_base.py +++ b/numpy/lib/twodim_base.py @@ -589,16 +589,18 @@ def histogram2d(x, y, bins=10, range=None, normed=False, weights=None): y : array_like, shape (N,) An array containing the y coordinates of the points to be histogrammed. - bins : int or [int, int] or array_like or [array, array], optional + bins : int or array_like or [int, int] or [array, array], optional The bin specification: * If int, the number of bins for the two dimensions (nx=ny=bins). - * If [int, int], the number of bins in each dimension - (nx, ny = bins). * If array_like, the bin edges for the two dimensions (x_edges=y_edges=bins). + * If [int, int], the number of bins in each dimension + (nx, ny = bins). * If [array, array], the bin edges in each dimension (x_edges, y_edges = bins). + * A combination [int, array] or [array, int], where int + is the number of bins and array is the bin edges. range : array_like, shape(2,2), optional The leftmost and rightmost edges of the bins along each dimension diff --git a/numpy/ma/core.py b/numpy/ma/core.py index 5d878363d..78f5fbc4d 100644 --- a/numpy/ma/core.py +++ b/numpy/ma/core.py @@ -3802,7 +3802,10 @@ class MaskedArray(ndarray): else: if m is not nomask: self._mask += m - ndarray.__iadd__(self._data, np.where(self._mask, 0, getdata(other))) + ndarray.__iadd__( + self._data, + np.where(self._mask, self.dtype.type(0), getdata(other)) + ) return self #.... def __isub__(self, other): @@ -3814,7 +3817,10 @@ class MaskedArray(ndarray): self._mask += m elif m is not nomask: self._mask += m - ndarray.__isub__(self._data, np.where(self._mask, 0, getdata(other))) + ndarray.__isub__( + self._data, + np.where(self._mask, self.dtype.type(0), getdata(other)) + ) return self #.... def __imul__(self, other): @@ -3826,7 +3832,10 @@ class MaskedArray(ndarray): self._mask += m elif m is not nomask: self._mask += m - ndarray.__imul__(self._data, np.where(self._mask, 1, getdata(other))) + ndarray.__imul__( + self._data, + np.where(self._mask, self.dtype.type(1), getdata(other)) + ) return self #.... def __idiv__(self, other): @@ -3841,7 +3850,10 @@ class MaskedArray(ndarray): other_data = np.where(dom_mask, fval, other_data) # self._mask = mask_or(self._mask, new_mask) self._mask |= new_mask - ndarray.__idiv__(self._data, np.where(self._mask, 1, other_data)) + ndarray.__idiv__( + self._data, + np.where(self._mask, self.dtype.type(1), other_data) + ) return self #.... def __ifloordiv__(self, other): @@ -3856,7 +3868,10 @@ class MaskedArray(ndarray): other_data = np.where(dom_mask, fval, other_data) # self._mask = mask_or(self._mask, new_mask) self._mask |= new_mask - ndarray.__ifloordiv__(self._data, np.where(self._mask, 1, other_data)) + ndarray.__ifloordiv__( + self._data, + np.where(self._mask, self.dtype.type(1), other_data) + ) return self #.... def __itruediv__(self, other): @@ -3871,7 +3886,10 @@ class MaskedArray(ndarray): other_data = np.where(dom_mask, fval, other_data) # self._mask = mask_or(self._mask, new_mask) self._mask |= new_mask - ndarray.__itruediv__(self._data, np.where(self._mask, 1, other_data)) + ndarray.__itruediv__( + self._data, + np.where(self._mask, self.dtype.type(1), other_data) + ) return self #... def __ipow__(self, other): @@ -3879,7 +3897,10 @@ class MaskedArray(ndarray): other_data = getdata(other) other_mask = getmask(other) with np.errstate(divide='ignore', invalid='ignore'): - ndarray.__ipow__(self._data, np.where(self._mask, 1, other_data)) + ndarray.__ipow__( + self._data, + np.where(self._mask, self.dtype.type(1), other_data) + ) invalid = np.logical_not(np.isfinite(self._data)) if invalid.any(): if self._mask is not nomask: diff --git a/numpy/ma/extras.py b/numpy/ma/extras.py index 9adb9d368..b6082180a 100644 --- a/numpy/ma/extras.py +++ b/numpy/ma/extras.py @@ -497,8 +497,9 @@ def average(a, axis=None, weights=None, returned=False): The average along the specified axis. When returned is `True`, return a tuple with the average as the first element and the sum of the weights as the second element. The return type is `np.float64` - if `a` is of integer type, otherwise it is of the same type as `a`. - If returned, `sum_of_weights` is of the same type as `average`. + if `a` is of integer type and floats smaller than `float64`, or the + input data-type, otherwise. If returned, `sum_of_weights` is always + `float64`. Examples -------- @@ -1920,12 +1921,12 @@ def polyfit(x, y, deg, rcond=None, full=False, w=None, cov=False): m = mask_or(m, getmask(w)) if m is not nomask: + not_m = ~m if w is not None: - w = ~m*w - else: - w = ~m - - return np.polyfit(x, y, deg, rcond, full, w, cov) + w = w[not_m] + return np.polyfit(x[not_m], y[not_m], deg, rcond, full, w, cov) + else: + return np.polyfit(x, y, deg, rcond, full, w, cov) polyfit.__doc__ = ma.doc_note(np.polyfit.__doc__, polyfit.__doc__) diff --git a/numpy/ma/tests/test_core.py b/numpy/ma/tests/test_core.py index 975b7aceb..1d4462306 100644 --- a/numpy/ma/tests/test_core.py +++ b/numpy/ma/tests/test_core.py @@ -1723,6 +1723,13 @@ class TestMaskedArrayInPlaceArithmetics(TestCase): xm[2] = masked self.intdata = (x, y, xm) self.floatdata = (x.astype(float), y.astype(float), xm.astype(float)) + self.othertypes = np.typecodes['AllInteger'] + np.typecodes['AllFloat'] + self.othertypes = [np.dtype(_).type for _ in self.othertypes] + self.uint8data = ( + x.astype(np.uint8), + y.astype(np.uint8), + xm.astype(np.uint8) + ) def test_inplace_addition_scalar(self): # Test of inplace additions @@ -1990,6 +1997,226 @@ class TestMaskedArrayInPlaceArithmetics(TestCase): assert_equal(a, [[1, 1], [3, 3]]) assert_equal(a.mask, [[0, 1], [0, 1]]) + def test_inplace_addition_scalar_type(self): + # Test of inplace additions + for t in self.othertypes: + with warnings.catch_warnings(record=True) as w: + warnings.filterwarnings("always") + (x, y, xm) = (_.astype(t) for _ in self.uint8data) + xm[2] = masked + x += t(1) + assert_equal(x, y + t(1)) + xm += t(1) + assert_equal(xm, y + t(1)) + + assert_equal(len(w), 0, "Failed on type=%s." % t) + + def test_inplace_addition_array_type(self): + # Test of inplace additions + for t in self.othertypes: + with warnings.catch_warnings(record=True) as w: + warnings.filterwarnings("always") + (x, y, xm) = (_.astype(t) for _ in self.uint8data) + m = xm.mask + a = arange(10, dtype=t) + a[-1] = masked + x += a + xm += a + assert_equal(x, y + a) + assert_equal(xm, y + a) + assert_equal(xm.mask, mask_or(m, a.mask)) + + assert_equal(len(w), 0, "Failed on type=%s." % t) + + def test_inplace_subtraction_scalar_type(self): + # Test of inplace subtractions + for t in self.othertypes: + with warnings.catch_warnings(record=True) as w: + warnings.filterwarnings("always") + (x, y, xm) = (_.astype(t) for _ in self.uint8data) + x -= t(1) + assert_equal(x, y - t(1)) + xm -= t(1) + assert_equal(xm, y - t(1)) + + assert_equal(len(w), 0, "Failed on type=%s." % t) + + + def test_inplace_subtraction_array_type(self): + # Test of inplace subtractions + for t in self.othertypes: + with warnings.catch_warnings(record=True) as w: + warnings.filterwarnings("always") + (x, y, xm) = (_.astype(t) for _ in self.uint8data) + m = xm.mask + a = arange(10, dtype=t) + a[-1] = masked + x -= a + xm -= a + assert_equal(x, y - a) + assert_equal(xm, y - a) + assert_equal(xm.mask, mask_or(m, a.mask)) + + assert_equal(len(w), 0, "Failed on type=%s." % t) + + def test_inplace_multiplication_scalar_type(self): + # Test of inplace multiplication + for t in self.othertypes: + with warnings.catch_warnings(record=True) as w: + warnings.filterwarnings("always") + (x, y, xm) = (_.astype(t) for _ in self.uint8data) + x *= t(2) + assert_equal(x, y * t(2)) + xm *= t(2) + assert_equal(xm, y * t(2)) + + assert_equal(len(w), 0, "Failed on type=%s." % t) + + def test_inplace_multiplication_array_type(self): + # Test of inplace multiplication + for t in self.othertypes: + with warnings.catch_warnings(record=True) as w: + warnings.filterwarnings("always") + (x, y, xm) = (_.astype(t) for _ in self.uint8data) + m = xm.mask + a = arange(10, dtype=t) + a[-1] = masked + x *= a + xm *= a + assert_equal(x, y * a) + assert_equal(xm, y * a) + assert_equal(xm.mask, mask_or(m, a.mask)) + + assert_equal(len(w), 0, "Failed on type=%s." % t) + + def test_inplace_floor_division_scalar_type(self): + # Test of inplace division + for t in self.othertypes: + with warnings.catch_warnings(record=True) as w: + warnings.filterwarnings("always") + (x, y, xm) = (_.astype(t) for _ in self.uint8data) + x = arange(10, dtype=t) * t(2) + xm = arange(10, dtype=t) * t(2) + xm[2] = masked + x //= t(2) + xm //= t(2) + assert_equal(x, y) + assert_equal(xm, y) + + assert_equal(len(w), 0, "Failed on type=%s." % t) + + def test_inplace_floor_division_array_type(self): + # Test of inplace division + for t in self.othertypes: + with warnings.catch_warnings(record=True) as w: + warnings.filterwarnings("always") + (x, y, xm) = (_.astype(t) for _ in self.uint8data) + m = xm.mask + a = arange(10, dtype=t) + a[-1] = masked + x //= a + xm //= a + assert_equal(x, y // a) + assert_equal(xm, y // a) + assert_equal( + xm.mask, + mask_or(mask_or(m, a.mask), (a == t(0))) + ) + + assert_equal(len(w), 0, "Failed on type=%s." % t) + + def test_inplace_division_scalar_type(self): + # Test of inplace division + for t in self.othertypes: + with warnings.catch_warnings(record=True) as w: + warnings.filterwarnings("always") + (x, y, xm) = (_.astype(t) for _ in self.uint8data) + x = arange(10, dtype=t) * t(2) + xm = arange(10, dtype=t) * t(2) + xm[2] = masked + + # May get a DeprecationWarning or a TypeError. + # + # This is a consequence of the fact that this is true divide + # and will require casting to float for calculation and + # casting back to the original type. This will only be raised + # with integers. Whether it is an error or warning is only + # dependent on how stringent the casting rules are. + # + # Will handle the same way. + try: + x /= t(2) + assert_equal(x, y) + except (DeprecationWarning, TypeError) as e: + warnings.warn(str(e)) + try: + xm /= t(2) + assert_equal(xm, y) + except (DeprecationWarning, TypeError) as e: + warnings.warn(str(e)) + + if issubclass(t, np.integer): + assert_equal(len(w), 2, "Failed on type=%s." % t) + else: + assert_equal(len(w), 0, "Failed on type=%s." % t) + + def test_inplace_division_array_type(self): + # Test of inplace division + for t in self.othertypes: + with warnings.catch_warnings(record=True) as w: + warnings.filterwarnings("always") + (x, y, xm) = (_.astype(t) for _ in self.uint8data) + m = xm.mask + a = arange(10, dtype=t) + a[-1] = masked + + # May get a DeprecationWarning or a TypeError. + # + # This is a consequence of the fact that this is true divide + # and will require casting to float for calculation and + # casting back to the original type. This will only be raised + # with integers. Whether it is an error or warning is only + # dependent on how stringent the casting rules are. + # + # Will handle the same way. + try: + x /= a + assert_equal(x, y / a) + except (DeprecationWarning, TypeError) as e: + warnings.warn(str(e)) + try: + xm /= a + assert_equal(xm, y / a) + assert_equal( + xm.mask, + mask_or(mask_or(m, a.mask), (a == t(0))) + ) + except (DeprecationWarning, TypeError) as e: + warnings.warn(str(e)) + + if issubclass(t, np.integer): + assert_equal(len(w), 2, "Failed on type=%s." % t) + else: + assert_equal(len(w), 0, "Failed on type=%s." % t) + + def test_inplace_pow_type(self): + # Test keeping data w/ (inplace) power + for t in self.othertypes: + with warnings.catch_warnings(record=True) as w: + warnings.filterwarnings("always") + # Test pow on scalar + x = array([1, 2, 3], mask=[0, 0, 1], dtype=t) + xx = x ** t(2) + xx_r = array([1, 2 ** 2, 3], mask=[0, 0, 1], dtype=t) + assert_equal(xx.data, xx_r.data) + assert_equal(xx.mask, xx_r.mask) + # Test ipow on scalar + x **= t(2) + assert_equal(x.data, xx_r.data) + assert_equal(x.mask, xx_r.mask) + + assert_equal(len(w), 0, "Failed on type=%s." % t) + #------------------------------------------------------------------------------ class TestMaskedArrayMethods(TestCase): diff --git a/numpy/ma/tests/test_extras.py b/numpy/ma/tests/test_extras.py index 8902d67dc..9e1e8ba38 100644 --- a/numpy/ma/tests/test_extras.py +++ b/numpy/ma/tests/test_extras.py @@ -734,6 +734,21 @@ class TestPolynomial(TestCase): for (a, a_) in zip((C, R, K, S, D), (c, r, k, s, d)): assert_almost_equal(a, a_) + def test_polyfit_with_masked_NaNs(self): + x = np.random.rand(10) + y = np.random.rand(20).reshape(-1, 2) + + x[0] = np.nan + y[-1,-1] = np.nan + x = x.view(MaskedArray) + y = y.view(MaskedArray) + x[0] = masked + y[-1,-1] = masked + + (C, R, K, S, D) = polyfit(x, y, 3, full=True) + (c, r, k, s, d) = np.polyfit(x[1:-1], y[1:-1,:], 3, full=True) + for (a, a_) in zip((C, R, K, S, D), (c, r, k, s, d)): + assert_almost_equal(a, a_) class TestArraySetOps(TestCase): diff --git a/numpy/random/mtrand/mtrand.pyx b/numpy/random/mtrand/mtrand.pyx index ed582767a..703e9ec28 100644 --- a/numpy/random/mtrand/mtrand.pyx +++ b/numpy/random/mtrand/mtrand.pyx @@ -3324,10 +3324,10 @@ cdef class RandomState: .. math:: P(x;scale) = \\frac{x}{scale^2}e^{\\frac{-x^2}{2 \\cdotp scale^2}} - The Rayleigh distribution arises if the wind speed and wind direction are - both gaussian variables, then the vector wind velocity forms a Rayleigh - distribution. The Rayleigh distribution is used to model the expected - output from wind turbines. + The Rayleigh distribution would arise, for example, if the East + and North components of the wind velocity had identical zero-mean + Gaussian distributions. Then the wind speed would have a Rayleigh + distribution. References ---------- diff --git a/numpy/testing/utils.py b/numpy/testing/utils.py index b9bdaeb87..1051288c2 100644 --- a/numpy/testing/utils.py +++ b/numpy/testing/utils.py @@ -1244,7 +1244,7 @@ def _assert_valid_refcount(op): assert_(sys.getrefcount(i) >= rc) -def assert_allclose(actual, desired, rtol=1e-7, atol=0, +def assert_allclose(actual, desired, rtol=1e-7, atol=0, equal_nan=False, err_msg='', verbose=True): """ Raises an AssertionError if two objects are not equal up to desired @@ -1266,6 +1266,8 @@ def assert_allclose(actual, desired, rtol=1e-7, atol=0, Relative tolerance. atol : float, optional Absolute tolerance. + equal_nan : bool, optional. + If True, NaNs will compare equal. err_msg : str, optional The error message to be printed in case of failure. verbose : bool, optional @@ -1289,7 +1291,8 @@ def assert_allclose(actual, desired, rtol=1e-7, atol=0, """ import numpy as np def compare(x, y): - return np.core.numeric._allclose_points(x, y, rtol=rtol, atol=atol) + return np.core.numeric.isclose(x, y, rtol=rtol, atol=atol, + equal_nan=equal_nan) actual, desired = np.asanyarray(actual), np.asanyarray(desired) header = 'Not equal to tolerance rtol=%g, atol=%g' % (rtol, atol) @@ -107,7 +107,7 @@ def get_version_info(): GIT_REVISION = "Unknown" if not ISRELEASED: - FULLVERSION += '.dev+' + GIT_REVISION[:7] + FULLVERSION += '.dev0+' + GIT_REVISION[:7] return FULLVERSION, GIT_REVISION diff --git a/site.cfg.example b/site.cfg.example index 4a59f10e2..cc92edb59 100644 --- a/site.cfg.example +++ b/site.cfg.example @@ -9,7 +9,7 @@ # The format of the file is that of the standard library's ConfigParser module. # -# http://www.python.org/doc/current/lib/module-ConfigParser.html +# http://docs.python.org/3/library/configparser.html # # Each section defines settings that apply to one particular dependency. Some of # the settings are general and apply to nearly any section and are defined here. @@ -118,6 +118,15 @@ # library_dirs = /opt/intel/mkl/10.0.1.014/lib/32/ # lapack_libs = mkl_lapack # mkl_libs = mkl, guide +# +# On win-64, the following options compiles numpy with the MKL library +# dynamically linked. +# [mkl] +# include_dirs = C:\Program Files (x86)\Intel\Composer XE 2015\mkl\include +# library_dirs = C:\Program Files (x86)\Intel\Composer XE 2015\mkl\lib\intel64 +# mkl_libs = mkl_core_dll, mkl_intel_lp64_dll, mkl_intel_thread_dll +# lapack_libs = mkl_lapack95_lp64 + # UMFPACK # ------- diff --git a/tools/swig/numpy.i b/tools/swig/numpy.i index 181625057..d9f1d48f1 100644 --- a/tools/swig/numpy.i +++ b/tools/swig/numpy.i @@ -3071,13 +3071,13 @@ * %numpy_typemaps(long double, NPY_LONGDOUBLE, int) */ -%#ifdef __cplusplus +#ifdef __cplusplus %include <std_complex.i> %numpy_typemaps(std::complex<float>, NPY_CFLOAT , int) %numpy_typemaps(std::complex<double>, NPY_CDOUBLE, int) -%#endif +#endif #endif /* SWIGPYTHON */ diff --git a/tools/swig/test/testArray.py b/tools/swig/test/testArray.py index 2ccca19b4..8d9c79772 100755 --- a/tools/swig/test/testArray.py +++ b/tools/swig/test/testArray.py @@ -382,4 +382,4 @@ if __name__ == "__main__": print("NumPy version", np.__version__) print() result = unittest.TextTestRunner(verbosity=2).run(suite) - sys.exit(len(result.errors) + len(result.failures)) + sys.exit(bool(result.errors + result.failures)) diff --git a/tools/swig/test/testFarray.py b/tools/swig/test/testFarray.py index 15683b70b..0037dc9b3 100755 --- a/tools/swig/test/testFarray.py +++ b/tools/swig/test/testFarray.py @@ -156,4 +156,4 @@ if __name__ == "__main__": print("NumPy version", np.__version__) print() result = unittest.TextTestRunner(verbosity=2).run(suite) - sys.exit(len(result.errors) + len(result.failures)) + sys.exit(bool(result.errors + result.failures)) diff --git a/tools/swig/test/testFortran.py b/tools/swig/test/testFortran.py index a42af950e..be4134d4e 100644 --- a/tools/swig/test/testFortran.py +++ b/tools/swig/test/testFortran.py @@ -170,4 +170,4 @@ if __name__ == "__main__": print("NumPy version", np.__version__) print() result = unittest.TextTestRunner(verbosity=2).run(suite) - sys.exit(len(result.errors) + len(result.failures)) + sys.exit(bool(result.errors + result.failures)) diff --git a/tools/swig/test/testMatrix.py b/tools/swig/test/testMatrix.py index af234e0e9..7127678f7 100755 --- a/tools/swig/test/testMatrix.py +++ b/tools/swig/test/testMatrix.py @@ -359,4 +359,4 @@ if __name__ == "__main__": print("NumPy version", np.__version__) print() result = unittest.TextTestRunner(verbosity=2).run(suite) - sys.exit(len(result.errors) + len(result.failures)) + sys.exit(bool(result.errors + result.failures)) diff --git a/tools/swig/test/testSuperTensor.py b/tools/swig/test/testSuperTensor.py index ff1f86df2..b7765ea0a 100644 --- a/tools/swig/test/testSuperTensor.py +++ b/tools/swig/test/testSuperTensor.py @@ -385,4 +385,4 @@ if __name__ == "__main__": print "NumPy version", np.__version__ print result = unittest.TextTestRunner(verbosity=2).run(suite) - sys.exit(len(result.errors) + len(result.failures)) + sys.exit(bool(result.errors + result.failures)) diff --git a/tools/swig/test/testTensor.py b/tools/swig/test/testTensor.py index a9390ebb1..61dc82090 100755 --- a/tools/swig/test/testTensor.py +++ b/tools/swig/test/testTensor.py @@ -399,4 +399,4 @@ if __name__ == "__main__": print("NumPy version", np.__version__) print() result = unittest.TextTestRunner(verbosity=2).run(suite) - sys.exit(len(result.errors) + len(result.failures)) + sys.exit(bool(result.errors + result.failures)) diff --git a/tools/swig/test/testVector.py b/tools/swig/test/testVector.py index e7d019cf7..eaaa75102 100755 --- a/tools/swig/test/testVector.py +++ b/tools/swig/test/testVector.py @@ -378,4 +378,4 @@ if __name__ == "__main__": print("NumPy version", np.__version__) print() result = unittest.TextTestRunner(verbosity=2).run(suite) - sys.exit(len(result.errors) + len(result.failures)) + sys.exit(bool(result.errors + result.failures)) |