Merge pull request #12485 from amiller27/remove-asserts

Replaced some asserts in non-test code with exceptions (#5090)

.ci/README.rst

@@ -0,0 +1,10 @@
+This folder currently holds ``durations.json`` and ``blacklisted.json`` which are used in our pytest configuration.
+``durations.json`` allows the two options ``--quickcheck`` and ``--veryquickcheck`` to be passed to ``pytest``.
+Currently, the list of slow tests is updated intermittently by running::
+
+  $ ./.ci/generate_durations_log.sh
+  $ ./.ci/parse_durations_log.py
+  $ git commit -am "pytest: Updated .ci/durations.json"
+
+make sure you have a C compiler, a Fortran compiler, numpy, scipy, cython, matplotlib & sage installed (the latter can
+be tricky, in which case just make sure that no sage tests are removed in the git diff).

.ci/generate_durations_log.sh

@@ -1,4 +1,4 @@
-#!/bin/bash
+#!/bin/bash -e
 ABS_REPO_PATH=$(unset CDPATH && cd "$(dirname "$0")/.." && echo $PWD)
 cat <<EOF >${ABS_REPO_PATH}/.ci/blacklisted.json
 {
@@ -10,4 +10,4 @@ cat <<EOF >${ABS_REPO_PATH}/.ci/blacklisted.json
     ]
 }
 EOF
-python3 -m pytest -ra --durations 0 --verbose | tee $ABS_REPO_PATH/.ci/durations.log
+${PYTHON:-python} -m pytest -ra --durations 0 --verbose | tee $ABS_REPO_PATH/.ci/durations.log

.github/FUNDING.yml

@@ -0,0 +1,3 @@
+# These are supported funding model platforms
+
+custom: ['https://numfocus.salsalabs.org/donate-to-sympy/index.html']

.gitignore

@@ -7,9 +7,12 @@
 # track a file that is ignored for some reason, you have to use
 # "git add -f". See "git help gitignore" for more information.
 
+# Virtualenv
+/.venv/
 
 # Regular Python bytecode file
 *.pyc
+__pycache__/
 
 # Optimized Python bytecode file
 *.pyo
@@ -17,6 +20,12 @@
 # Vim's swap files
 *.sw[op]
 
+# VSCode settings file
+.vscode/
+
+# PyCharm settings file
+/.idea/
+
 # Generated files from Jython
 *$py.class
 

.travis.yml

@@ -8,7 +8,7 @@ env:
   matrix:
   - SPLIT="1/2" TEST_SYMPY="true"
   - SPLIT="2/2" TEST_SYMPY="true"
-  - TEST_DOCTESTS="true" FASTCACHE="false" TEST_SETUP="true"
+  - TEST_DOCTESTS="true" FASTCACHE="false" TEST_SETUP="true" TEST_EXAMPLES="true"
   global:
     - secure: "YIEZal9EBTL+fg2YmoZoS8Bvt3eAVUOZjb38CtqpzR2CCSXWoUk35KG23m2NknlY1iKfYJyt7XWBszT/VKOQEbWQq7PIakV4vIByrWacgBxy1x3WC+rZoW7TX+JJiL+y942qIYbMoNMMB8xFpE5RDLSjSecMpFhJJXoafVTvju8="
 dist: trusty
@@ -31,7 +31,7 @@ matrix:
       dist: xenial
       sudo: true
       env:
-        - TEST_DOCTESTS="true" FASTCACHE="false" TEST_SETUP="true"
+        - TEST_DOCTESTS="true" FASTCACHE="false" TEST_SETUP="true" TEST_EXAMPLES="true"
     - python: 3.7
       dist: xenial
       sudo: true
@@ -44,11 +44,26 @@ matrix:
         - SPLIT="2/2" TEST_SYMPY="true"
 
     - stage: test
-      python: 2.7
+      python: 3.8-dev
+      dist: xenial
+      sudo: true
+      env:
+        - TEST_DOCTESTS="true" FASTCACHE="false" TEST_SETUP="true" TEST_EXAMPLES="true"
+    - python: 3.8-dev
+      dist: xenial
+      sudo: true
+      env:
+        - SPLIT="1/2" TEST_SYMPY="true"
+    - python: 3.8-dev
+      dist: xenial
+      sudo: true
+      env:
+        - SPLIT="2/2" TEST_SYMPY="true"
+    - python: 2.7
       env:
         - TEST_ASCII="true"
         # space separated list of optional dependencies(conda packages) to install and test
-        - TEST_OPT_DEPENDENCY="numpy scipy gmpy2 matplotlib>=2.2 theano llvmlite autowrap cython wurlitzer python-symengine=0.3.* tensorflow numexpr ipython antlr-python-runtime>=4.7,<4.8 antlr>=4.7,<4.8 cloudpickle python=2.7 pyglet"
+        - TEST_OPT_DEPENDENCY="numpy scipy gmpy2 matplotlib>=2.2 theano llvmlite autowrap cython wurlitzer python-symengine=0.3.* tensorflow numexpr ipython antlr-python-runtime>=4.7,<4.8 antlr>=4.7,<4.8 cloudpickle python=2.7 pyglet pycosat python-clang"
       addons:
         apt:
           packages:
@@ -63,7 +78,7 @@ matrix:
       # and we aren't actually using the Travis Python anyway.
       env:
         - TEST_ASCII="true"
-        - TEST_OPT_DEPENDENCY="numpy scipy gmpy2 matplotlib theano llvmlite autowrap cython wurlitzer python-symengine=0.3.* tensorflow numexpr ipython antlr-python-runtime>=4.7,<4.8 antlr>=4.7,<4.8 cloudpickle pyglet"
+        - TEST_OPT_DEPENDENCY="numpy scipy gmpy2 matplotlib theano llvmlite autowrap cython wurlitzer python-symengine=0.3.* tensorflow numexpr ipython antlr-python-runtime>=4.7,<4.8 antlr>=4.7,<4.8 cloudpickle pyglet pycosat lfortran python-clang"
         - TEST_SAGE="true"
       addons:
         apt:
@@ -151,7 +166,7 @@ matrix:
     # duplicated here and in that section below.
     - python: "pypy"
       env:
-        - TEST_DOCTESTS="true"
+        - TEST_DOCTESTS="true" TEST_EXAMPLES="true"
         - FASTCACHE="false"
       addons:
         apt:
@@ -229,6 +244,9 @@ before_install:
       pip install coverage;
       coverage debug sys;
     fi
+  - if [ "$TRAVIS_PYTHON_VERSION" == "3.8-dev" ]; then
+        pip uninstall -y Cython;
+    fi
   - |
     if [[ -n "${TEST_OPT_DEPENDENCY}" ]]; then
     # We do this conditionally because it saves us some downloading if the
@@ -260,14 +278,13 @@ before_install:
             conda clean --all;
         fi
         if [[ "$TRAVIS_PYTHON_VERSION" == "3.6" ]]; then
-            # Change this back to conda once matchpy 0.5.0 gets updated on conda-forge
-            pip install matchpy;
+            conda install matchpy;
         fi
     elif [ "$TRAVIS_PYTHON_VERSION" != "pypy" ]; then
         pip list | grep "numpy" && pip uninstall -y numpy;
     fi
   - if [[ "${TEST_SPHINX}" == "true" ]]; then
-      pip install "sphinx" "docutils" doctr matplotlib;
+      pip install "sphinx" "docutils" doctr matplotlib sphinx-math-dollar;
     fi
 
 install:
@@ -282,6 +299,9 @@ install:
   - pip install mpmath
   # -We:invalid makes invalid escape sequences error in Python 3.6. See
   # -#12028.
+  # SyntaxWarning flag for catching errors in Python3.8
+  # Issue -  #16973. -We:invalid can be dropped from 3.8 onwards, but
+  # it needs to be there for earlier versions.
   - |
     if [[ "${TEST_SETUP}" == "true" ]]; then
       # The install cycle below is to test installation on systems without
@@ -291,7 +311,7 @@ install:
       ~/.venv-no-setuptools/bin/pip uninstall -y setuptools;
       ~/.venv-no-setuptools/bin/python -We:invalid setup.py -q install;
     fi
-    python -We:invalid -m compileall -f -q sympy/;
+    python -We:invalid -We::SyntaxWarning -m compileall -f -q sympy/;
     python -We:invalid setup.py -q install;
     pip list --format=columns;
 

bin/test_external_imports.py

@@ -25,7 +25,7 @@ hard_dependencies = ['mpmath']
 # These libraries are optional, but are always imported at SymPy import time
 # when they are installed. External libraries should only be added to this
 # list if they are required for core SymPy functionality.
-hard_optional_dependencies = ['gmpy', 'gmpy2', 'fastcache']
+hard_optional_dependencies = ['gmpy', 'gmpy2', 'fastcache', 'pycosat']
 
 import sys
 import os

bin/test_travis.sh

@@ -137,9 +137,8 @@ test_list = [
     # ipython
     '*ipython*',
 
-    # antlr
-    'sympy/parsing/tests/test_autolev',
-    'sympy/parsing/tests/test_latex',
+    # antlr, lfortran, clang
+    'sympy/parsing/',
 
     # matchpy
     '*rubi*',
@@ -151,6 +150,11 @@ test_list = [
 
     # cloudpickle
     'pickling',
+
+    # pycosat
+    'sympy/logic',
+    'sympy/assumptions',
+
 ]
 
 blacklist = [
@@ -180,15 +184,19 @@ doctest_list = [
     # ipython
     '*ipython*',
 
-    # antlr
-    'sympy/parsing/autolev',
-    'sympy/parsing/latex',
+    # antlr, lfortran, clang
+    'sympy/parsing/',
 
     # matchpy
     '*rubi*',
 
     # codegen
     'sympy/codegen/',
+
+    # pycosat
+    'sympy/logic',
+    'sympy/assumptions',
+
 ]
 
 if not (sympy.test(*test_list, blacklist=blacklist) and sympy.doctest(*doctest_list)):
@@ -242,3 +250,8 @@ fi
 if [[ "${TEST_COVERAGE}" == "true" ]]; then
     unset COVERAGE_PROCESS_START
 fi
+
+if [[ "${TEST_EXAMPLES}" == "true" ]]; then
+    # No need to change directory if executed after the rst doctest
+    examples/all.py -q
+fi

data/TeXmacs/bin/tm_sympy

@@ -11,7 +11,7 @@
 # heuristics to know where a multi-line expression is
 # broken. As a result you can't use more than one space
 # in a sequence. However you can and must indent your
-# expression using standard Pyhon rules.
+# expression using standard Python rules.
 #
 # You can retrieve the last output using '_' built-in
 # symbol. If the previous command did not generate

doc/README.rst

@@ -7,6 +7,7 @@ Debian/Ubuntu
 To make the html documentation, install the prerequisites::
 
     apt-get install python-sphinx texlive-latex-recommended dvipng librsvg2-bin imagemagick docbook2x graphviz
+    python -m pip install sphinx-math-dollar
 
 and do::
 
@@ -30,6 +31,7 @@ Fedora
 Fedora (and maybe other RPM based distributions), install the prerequisites::
 
     dnf install python3-sphinx librsvg2 ImageMagick docbook2X texlive-dvipng-bin texlive-scheme-medium librsvg2-tools
+    python -m pip install sphinx-math-dollar
 
 After that, run::
 
@@ -47,6 +49,30 @@ And view it at::
 
     _build/html/index.html
 
+Mac
+---
+
+- Install MacTex: http://www.tug.org/mactex/
+- Install homebrew: https://brew.sh/
+
+Install these packages with homebrew:
+
+```bash
+brew install imagemagick graphviz docbook librsvg
+```
+
+Install these packages with either pip or conda:
+
+```bash
+python -m pip install mpmath matplotlib sphinx sphinx-math-dollar
+```
+
+or
+
+```bash
+conda install -c conda-forge mpmpath matplotlib sphinx sphinx-math-dollar
+```
+
 Windows 10
 ----------
 

doc/src/_static/default.css_t

@@ -9,7 +9,7 @@
  *
  */
 
-@import url("basic.css");
+@import url("classic.css");
 @import url(https://fonts.googleapis.com/css?family=Open+Sans);
 
 /* -- page layout ----------------------------------------------------------- */

doc/src/conf.py

@@ -12,6 +12,9 @@
 # serve to show the default value.
 
 import sys
+import inspect
+import os
+import subprocess
 import sympy
 
 # If your extensions are in another directory, add it here.
@@ -22,7 +25,7 @@ sys.path = ['ext'] + sys.path
 
 # Add any Sphinx extension module names here, as strings. They can be extensions
 # coming with Sphinx (named 'sphinx.addons.*') or your custom ones.
-extensions = ['sphinx.ext.autodoc', 'sphinx.ext.viewcode',
+extensions = ['sphinx.ext.autodoc', 'sphinx.ext.linkcode', 'sphinx_math_dollar',
               'sphinx.ext.mathjax', 'numpydoc', 'sympylive',
               'sphinx.ext.graphviz', 'matplotlib.sphinxext.plot_directive']
 
@@ -204,5 +207,83 @@ texinfo_documents = [
 ]
 
 # Use svg for graphviz
-
 graphviz_output_format = 'svg'
+
+
+# Requried for linkcode extension.
+# Get commit hash from the external file.
+commit_hash_filepath = '../commit_hash.txt'
+commit_hash = None
+if os.path.isfile(commit_hash_filepath):
+    with open(commit_hash_filepath, 'r') as f:
+        commit_hash = f.readline()
+
+# Get commit hash from the external file.
+if not commit_hash:
+    try:
+        commit_hash = subprocess.check_output(['git', 'rev-parse', 'HEAD'])
+        commit_hash = commit_hash.decode('ascii')
+        commit_hash = commit_hash.rstrip()
+    except:
+        import warnings
+        warnings.warn(
+            "Failed to get the git commit hash as the command " \
+            "'git rev-parse HEAD' is not working. The commit hash will be " \
+            "assumed as the SymPy master, but the lines may be misleading " \
+            "or nonexistent as it is not the correct branch the doc is " \
+            "built with. Check your installation of 'git' if you want to " \
+            "resolve this warning.")
+        commit_hash = 'master'
+
+fork = 'sympy'
+blobpath = \
+    "https://github.com/{}/sympy/blob/{}/sympy/".format(fork, commit_hash)
+
+
+def linkcode_resolve(domain, info):
+    """Determine the URL corresponding to Python object."""
+    if domain != 'py':
+        return
+
+    modname = info['module']
+    fullname = info['fullname']
+
+    submod = sys.modules.get(modname)
+    if submod is None:
+        return
+
+    obj = submod
+    for part in fullname.split('.'):
+        try:
+            obj = getattr(obj, part)
+        except Exception:
+            return
+
+    # strip decorators, which would resolve to the source of the decorator
+    # possibly an upstream bug in getsourcefile, bpo-1764286
+    try:
+        unwrap = inspect.unwrap
+    except AttributeError:
+        pass
+    else:
+        obj = unwrap(obj)
+
+    try:
+        fn = inspect.getsourcefile(obj)
+    except Exception:
+        fn = None
+    if not fn:
+        return
+
+    try:
+        source, lineno = inspect.getsourcelines(obj)
+    except Exception:
+        lineno = None
+
+    if lineno:
+        linespec = "#L%d-L%d" % (lineno, lineno + len(source) - 1)
+    else:
+        linespec = ""
+
+    fn = os.path.relpath(fn, start=os.path.dirname(sympy.__file__))
+    return blobpath + fn + linespec

doc/src/guide.rst

@@ -127,7 +127,7 @@ For easier use, there is a syntactic sugar for expressions.
 For example, ``cos(x) + 1`` is equal to ``cos(x).__add__(1)`` is equal to
 ``Add(cos(x), Integer(1))``
 
-This works becauses Python allows you to override the behavior of built-in
+This works because Python allows you to override the behavior of built-in
 operations ``+, *, -, /`` etc., by defining methods with names ``__add__``,
 ``__mul__`` etc., so ``cos(x) + 1`` is evaluated as ``cos(x).__add__(1)``.
 

doc/src/modules/algebras.rst

@@ -1,6 +1,6 @@
-===============
-Algebras Module
-===============
+========
+Algebras
+========
 
 Introduction
 ------------

doc/src/modules/assumptions/handlers/common.rst

@@ -0,0 +1,6 @@
+======
+Common
+======
+
+.. automodule:: sympy.assumptions.handlers.common
+   :members:

doc/src/modules/assumptions/handlers/index.rst

@@ -11,6 +11,8 @@ Contents
     :maxdepth: 3
 
     calculus.rst
+    common.rst
+    matrices.rst
     ntheory.rst
     order.rst
     sets.rst

doc/src/modules/assumptions/handlers/matrices.rst

@@ -0,0 +1,6 @@
+========
+Matrices
+========
+
+.. automodule:: sympy.assumptions.handlers.matrices
+   :members:

doc/src/modules/assumptions/index.rst

@@ -1,8 +1,8 @@
 .. _assumptions_module:
 
-==================
-Assumptions Module
-==================
+===========
+Assumptions
+===========
 
 .. automodule:: sympy.assumptions
 

doc/src/modules/calculus/index.rst

@@ -10,6 +10,8 @@ Calculus
 .. automodule:: sympy.calculus.singularities
     :members:
 
+.. _finite_diff:
+
 .. automodule:: sympy.calculus.finite_diff
     :members:
 

doc/src/modules/categories.rst

@@ -1,5 +1,6 @@
-Category Theory Module
-======================
+===============
+Category Theory
+===============
 
 Introduction
 ------------

doc/src/modules/codegen.rst

@@ -1,8 +1,8 @@
 .. _codegen_prose:
 
-================================================
-Structural Details of Code Generation with SymPy
-================================================
+===============
+Code Generation
+===============
 
 Several submodules in SymPy allow one to generate directly compilable and
 executable code in a variety of different programming languages from SymPy
@@ -164,8 +164,8 @@ An example of Mathematica code printer::
 
     >>> expr = summation(expr, (n, -1, 1))
     >>> mathematica_code(expr)
-    T*x[-T]*Sin[(T + t)/T]/(T + t) + T*x[T]*Sin[(-T + t)/T]/(-T + t) + T*x[0]*Sin[
-    t/T]/t
+    T*(x[-T]*Sin[(T + t)/T]/(T + t) + x[T]*Sin[(-T + t)/T]/(-T + t) + x[0]*Sin[t/T
+    ]/t)
 
 We can go through a common expression in different languages we support and see
 how it works::
@@ -548,6 +548,9 @@ Classes and functions for rewriting expressions (sympy.codegen.rewriting)
 .. automodule:: sympy.codegen.rewriting
    :members:
 
+.. automodule:: sympy.codegen.matrix_nodes
+   :members:
+
 Tools for simplifying expressions using approximations (sympy.codegen.approximations)
 -------------------------------------------------------------------------------------
 
@@ -612,3 +615,10 @@ Fortran utilities (sympy.codegen.futils)
 
 .. automodule:: sympy.codegen.futils
    :members:
+
+
+Array utilities (sympy.codegen.array_utils)
+-------------------------------------------
+
+.. automodule:: sympy.codegen.array_utils
+   :members:

doc/src/modules/combinatorics/index.rst

@@ -1,8 +1,8 @@
 .. _combinatiorics-docs:
 
-====================
-Combinatorics Module
-====================
+=============
+Combinatorics
+=============
 
 Contents
 ========
@@ -23,3 +23,4 @@ Contents
    testutil.rst
    tensor_can.rst
    fp_groups.rst
+   pc_groups.rst

doc/src/modules/combinatorics/pc_groups.rst

@@ -0,0 +1,307 @@
+Polycyclic Groups
+=================
+
+Introduction
+------------
+
+This module presents the functionality designed for computing with
+polycyclic groups(PcGroup for short). The name of the corresponding
+SymPy object is ``PolycyclicGroup``. The functions or classes described
+here are studied under **Computational Group Theory**.
+
+Overview of functionalities
+```````````````````````````
+
+* The construction of PolycyclicGroup from a given PermutationGroup.
+
+* Computation of polycyclic generating sequence(pcgs for short) and
+  polycyclic series(pc_series).
+
+* Computation of relative order for polycyclic series.
+
+* Implementation of class Collector which can be treated as a base for polycylic groups.
+
+* Implementation of polycyclic group presentation(pc_presentation for short).
+
+* Computation of exponent vector, depth and leading exponent for a given element
+  of a polycyclic group.
+
+For a description of fundamental algorithms of polycyclic groups, we
+often make use of *Handbook of Computational Group Theory*.
+
+
+The Construction of Polycyclic Groups
+-------------------------------------
+
+Given a Permutation Group, A Polycyclic Group is constructed by computing the
+corresponding polycylic generating sequence, polycyclic series and it's
+relative order.
+
+Attributes of PolycyclicGroup
+`````````````````````````````
+
+* ``pc_sequence`` : Polycyclic sequence is formed by collecting all the missing 
+  generators between the adjacent groups in the derived series of given
+  permutation group.
+
+* ``pc_series`` : Polycyclic series is formed by adding all the missing generators
+  of ``der[i+1]`` in ``der[i]``, where ``der`` represents derived series.
+
+* ``relative_order`` : A list, computed by the ratio of adjacent groups in pc_series.
+
+* ``collector`` : By default, it is None. Collector class provides the polycyclic presentation.
+
+>>> from sympy.combinatorics.named_groups import SymmetricGroup
+>>> G = SymmetricGroup(4)
+>>> PcGroup = G.polycyclic_group()
+>>> len(PcGroup.pcgs)
+4
+>>> pc_series = PcGroup.pc_series
+>>> pc_series[0]==G
+True
+>>> gen = pc_series[len(pc_series)-1].generators[0]
+>>> gen.is_identity
+True
+>>> PcGroup.relative_order
+[2, 3, 2, 2]
+
+
+The Construction of Collector
+-----------------------------
+
+Collector is one of the attributes of class PolycyclicGroup.
+
+Attributes of Collector
+```````````````````````
+
+Collector posses all the attributes of PolycyclicGroup, In addition there are
+few more attributes which are defined below:
+
+* ``free_group`` : free_group provides the mapping of polycyclic generating sequence with
+  the free group elements.
+
+* ``pc_presentation`` : Provides the presentation of polycyclic groups with the
+  help of power and conjugate relators.
+
+>>> from sympy.combinatorics.named_groups import SymmetricGroup
+>>> G = SymmetricGroup(3)
+>>> PcGroup = G.polycyclic_group()
+>>> Collector = PcGroup.collector
+>>> Collector.free_group
+<free group on the generators (x0, x1)>
+>>> Collector.pc_presentation
+{x0**2: (), x1**3: (), x0**-1*x1*x0: x1**2}
+
+
+Computation of Minimal Uncollected Subword
+``````````````````````````````````````````
+
+A word ``V`` defined on generators in the free_group of pc_group is a minimal
+uncollected subword of the word ``W`` if ``V`` is a subword of ``W`` and it has one of
+the following form:
+
+* `v = {x_{i+1}}^{a_j}x_i`
+* `v = {x_{i+1}}^{a_j}{x_i}^{-1}`
+* `v = {x_i}^{a_j}`
+
+`a_j \notin \{0, \ldots \mathrm{relative\_order}[j]-1\}`.
+
+>>> from sympy.combinatorics.named_groups import SymmetricGroup
+>>> from sympy.combinatorics.free_groups import free_group
+>>> G = SymmetricGroup(4)
+>>> PcGroup = G.polycyclic_group()
+>>> collector = PcGroup.collector
+>>> F, x1, x2 = free_group("x1, x2")
+>>> word = x2**2*x1**7
+>>> collector.minimal_uncollected_subword(word)
+((x2, 2),)
+
+
+Computation of Subword Index
+````````````````````````````
+
+For a given word and it's subword, subword_index computes the
+starting and ending index of the subword in the word.
+
+>>> from sympy.combinatorics.named_groups import SymmetricGroup
+>>> from sympy.combinatorics.free_groups import free_group
+>>> G = SymmetricGroup(4)
+>>> PcGroup = G.polycyclic_group()
+>>> collector = PcGroup.collector
+>>> F, x1, x2 = free_group("x1, x2")
+>>> word = x2**2*x1**7
+>>> w = x2**2*x1
+>>> collector.subword_index(word, w)
+(0, 3)
+>>> w = x1**7
+>>> collector.subword_index(word, w)
+(2, 9)
+
+
+Computation of Collected Word
+`````````````````````````````
+
+A word ``W`` is called collected, if ``W`` `= {x_{i_1}}^{a_1} \ldots {x_{i_r}}^{a_r}`
+with `i_1 < i_2< \ldots < i_r` and `a_j` is in `\{1 \ldots s_{j-1}\}`, where `s_j`
+represents the respective relative order.
+
+>>> from sympy.combinatorics.named_groups import SymmetricGroup
+>>> from sympy.combinatorics.perm_groups import PermutationGroup
+>>> from sympy.combinatorics.free_groups import free_group
+>>> G = SymmetricGroup(4)
+>>> PcGroup = G.polycyclic_group()
+>>> collector = PcGroup.collector
+>>> F, x0, x1, x2, x3 = free_group("x0, x1, x2, x3")
+>>> word = x3*x2*x1*x0
+>>> collected_word = collector.collected_word(word)
+>>> free_to_perm = {}
+>>> free_group = collector.free_group
+>>> for sym, gen in zip(free_group.symbols, collector.pcgs):
+...     free_to_perm[sym] = gen
+>>> G1 = PermutationGroup()
+>>> for w in word:
+...     sym = w[0]
+...     perm = free_to_perm[sym]
+...     G1 = PermutationGroup([perm] + G1.generators)
+>>> G2 = PermutationGroup()
+>>> for w in collected_word:
+...     sym = w[0]
+...     perm = free_to_perm[sym]
+...     G2 = PermutationGroup([perm] + G2.generators)
+>>> G1 == G2
+True
+
+
+Computation of Polycyclic Presentation
+--------------------------------------
+
+The computation of presentation starts from the bottom of the pcgs and polycyclic series.
+Storing all the previous generators from pcgs and then taking the last generator
+as the generator which acts as a conjugator and conjugates all the previous
+generators in the list.
+
+To get a clear picture, start with an example of SymmetricGroup(4). For S(4) there are 4
+generators in pcgs say `[x_0, x_1, x_2, x_3]` and the relative_order vector is [2, 3, 2, 2].
+Starting from bottom of this sequence the presentation is computed in order as below.
+
+using only `[x_3]` from ``pcgs`` and ``pc_series[4]`` compute:
+
+* `x_3^2`
+
+using only `[x_3]` from ``pcgs`` and ``pc_series[3]`` compute:
+
+* `x_2^2`
+* `x_2^{-1}x_3x_2`
+
+using `[x_3, x_2]` from ``pcgs`` and ``pc_series[2]`` compute:
+
+* `x_1^3`
+* `x_1^{-1}x_3x_1`
+* `x_1^{-1}x_2x_1`
+
+using `[x_3, x_2, x_1]` from ``pcgs`` and ``pc_series[1]`` compute:
+
+* `x_0^2`
+* `x_0^{-1}x_3x_0`
+* `x_0^{-1}x_2x_0`
+* `x_0^{-1}x_1x_0`
+
+
+One thing to note is same group can have different pcgs due to variying derived_series which,
+results in different polycyclic presentations.
+
+>>> from sympy.combinatorics.named_groups import SymmetricGroup
+>>> from sympy.combinatorics.permutations import Permutation
+>>> G = SymmetricGroup(4)
+>>> PcGroup = G.polycyclic_group()
+>>> collector = PcGroup.collector
+>>> pcgs = PcGroup.pcgs
+>>> len(pcgs)
+4
+>>> free_group = collector.free_group
+>>> pc_resentation = collector.pc_presentation
+>>> free_to_perm = {}
+>>> for s, g in zip(free_group.symbols, pcgs):
+...     free_to_perm[s] = g
+>>> for k, v in pc_resentation.items():
+...     k_array = k.array_form
+...     if v != ():
+...        v_array = v.array_form
+...     lhs = Permutation()
+...     for gen in k_array:
+...         s = gen[0]
+...         e = gen[1]
+...         lhs = lhs*free_to_perm[s]**e
+...     if v == ():
+...         assert lhs.is_identity
+...         continue
+...     rhs = Permutation()
+...     for gen in v_array:
+...         s = gen[0]
+...         e = gen[1]
+...         rhs = rhs*free_to_perm[s]**e
+...     assert lhs == rhs
+
+
+Computation of Exponent Vector
+------------------------------
+
+Any generator of the polycyclic group can be represented with the help of it's
+polycyclic generating sequence. Hence, the length of exponent vector is equal to
+the length of the pcgs.
+
+A given generator ``g`` of the polycyclic group, can be represented as
+`g = x_1^{e_1} \ldots x_n^{e_n}`, where `x_i` represents polycyclic generators 
+and ``n`` is the number of generators in the free_group equal to the length of pcgs.
+
+>>> from sympy.combinatorics.named_groups import SymmetricGroup
+>>> from sympy.combinatorics.permutations import Permutation
+>>> G = SymmetricGroup(4)
+>>> PcGroup = G.polycyclic_group()
+>>> collector = PcGroup.collector
+>>> pcgs = PcGroup.pcgs
+>>> collector.exponent_vector(G[0])
+[1, 0, 0, 0]
+>>> exp = collector.exponent_vector(G[1])
+>>> g = Permutation()
+>>> for i in range(len(exp)):
+...     g = g*pcgs[i]**exp[i] if exp[i] else g
+>>> assert g == G[1]
+
+
+Depth of Polycyclic generator
+`````````````````````````````
+
+Depth of a given polycyclic generator is defined as the index of the first
+non-zero entry in the exponent vector.
+
+>>> from sympy.combinatorics.named_groups import SymmetricGroup
+>>> G = SymmetricGroup(3)
+>>> PcGroup = G.polycyclic_group()
+>>> collector = PcGroup.collector
+>>> collector.depth(G[0])
+2
+>>> collector.depth(G[1])
+1
+
+
+Computation of Leading Exponent
+```````````````````````````````
+
+Leading exponent represents the exponent of polycyclic generator at the above depth.
+
+>>> from sympy.combinatorics.named_groups import SymmetricGroup
+>>> G = SymmetricGroup(3)
+>>> PcGroup = G.polycyclic_group()
+>>> collector = PcGroup.collector
+>>> collector.leading_exponent(G[1])
+1
+
+
+Bibliography
+------------
+
+.. [Ho05] Derek F. Holt,
+    Handbook of Computational Group Theory.
+    In the series 'Discrete Mathematics and its Applications',
+    `Chapman & Hall/CRC 2005, xvi + 514 p <https://www.crcpress.com/Handbook-of-Computational-Group-Theory/Holt-Eick-OBrien/p/book/9781584883722>`_.

doc/src/modules/concrete.rst

@@ -1,5 +1,6 @@
-Concrete Mathematics
-====================
+========
+Concrete
+========
 
 Hypergeometric terms
 --------------------

doc/src/modules/core.rst

@@ -1,5 +1,6 @@
-SymPy Core
-==========
+====
+Core
+====
 
 sympify
 -------

doc/src/modules/crypto.rst

@@ -1,5 +1,6 @@
-Basic Cryptography Module
-=========================
+============
+Cryptography
+============
 
 .. warning::
 
@@ -141,3 +142,7 @@ substitutions at different times in the message.
 .. autofunction:: encipher_gm
 
 .. autofunction:: decipher_gm
+
+.. autofunction:: encipher_railfence
+
+.. autofunction:: decipher_railfence

doc/src/modules/diffgeom.rst

@@ -1,5 +1,6 @@
-Differential Geometry Module
-============================
+=====================
+Differential Geometry
+=====================
 
 .. module:: sympy.diffgeom
 

doc/src/modules/discrete.rst

@@ -1,5 +1,6 @@
-Discrete Module
-===============
+========
+Discrete
+========
 
 The ``discrete`` module in SymPy implements methods to compute discrete
 transforms and convolutions of finite sequences.

doc/src/modules/evalf.rst

@@ -1,6 +1,7 @@
 .. _evalf-label:
 
-Numerical evaluation
+====================
+Numerical Evaluation
 ====================
 
 Basics

doc/src/modules/functions/index.rst

@@ -1,7 +1,8 @@
 .. _functions:
 
-Functions Module
-****************
+=========
+Functions
+=========
 
 .. module:: sympy.functions
 

doc/src/modules/geometry/index.rst

@@ -1,6 +1,6 @@
-===============
-Geometry Module
-===============
+========
+Geometry
+========
 
 
 Introduction

doc/src/modules/holonomic/about.rst

@@ -36,7 +36,7 @@ Here `D` is the Differential Operator and `L` is called the annihilator
 of the function.
 
 A unique holonomic function can be defined from the annihilator and a set of
-initial conditons.
+initial conditions.
 For instance:
 
 .. math::

doc/src/modules/holonomic/convert.rst

@@ -1,7 +1,7 @@
 Converting other representations to holonomic
 =============================================
 
-.. module:: sympy.holonomic.holonomic
+.. currentmodule:: sympy.holonomic.holonomic
 
 Converting hypergeometric functions
 -----------------------------------

doc/src/modules/holonomic/index.rst

@@ -1,8 +1,8 @@
 .. _holonomic-docs:
 
-===================
-Holonomic Functions
-===================
+=========
+Holonomic
+=========
 
 .. automodule:: sympy.holonomic
 

doc/src/modules/holonomic/operations.rst

@@ -29,7 +29,7 @@ a holonomic functions.
     >>> p * q
     HolonomicFunction((2) + (-2)*Dx + (1)*Dx**2, x, 0, [0, 1])
 
-.. module:: sympy.holonomic.holonomic
+.. currentmodule:: sympy.holonomic.holonomic
 
 Integration and Differentiation
 -------------------------------

doc/src/modules/holonomic/represent.rst

@@ -1,7 +1,7 @@
 Representation of holonomic functions in SymPy
 ==============================================
 
-.. module:: sympy.holonomic.holonomic
+.. currentmodule:: sympy.holonomic.holonomic
 
 Class :class:`DifferentialOperator` is used to represent the annihilator
 but we create differential operators easily using the function

doc/src/modules/holonomic/uses.rst

@@ -1,4 +1,4 @@
-.. module:: sympy.holonomic.holonomic
+.. currentmodule:: sympy.holonomic.holonomic
 
 Uses and Current limitations
 ============================

doc/src/modules/index.rst

@@ -16,49 +16,49 @@ access any SymPy module, or use the list below:
 .. toctree::
    :maxdepth: 2
 
-   core.rst
-   combinatorics/index.rst
-   ntheory.rst
-   crypto.rst
-   concrete.rst
-   discrete.rst
-   evalf.rst
+   algebras.rst
+   assumptions/index.rst
+   calculus/index.rst
+   categories.rst
    codegen.rst
-   numeric-computation.rst
+   combinatorics/index.rst
+   concrete.rst
+   core.rst
+   crypto.rst
+   diffgeom.rst
+   solvers/diophantine.rst
+   discrete.rst
    functions/index.rst
    geometry/index.rst
    holonomic/index.rst
+   simplify/hyperexpand.rst
+   solvers/inequalities.rst
    integrals/integrals.rst
    liealgebras/index.rst
    logic.rst
    matrices/index.rst
+   ntheory.rst
+   numeric-computation.rst
+   evalf.rst
+   solvers/ode.rst
+   parsing.rst
+   solvers/pde.rst
+   physics/index.rst
+   plotting.rst
    polys/index.rst
    printing.rst
-   plotting.rst
-   assumptions/index.rst
-   rewriting.rst
    series/index.rst
    sets.rst
    simplify/simplify.rst
-   simplify/hyperexpand.rst
-   stats.rst
-   solvers/ode.rst
-   solvers/pde.rst
    solvers/solvers.rst
-   solvers/diophantine.rst
-   solvers/inequalities.rst
    solvers/solveset.rst
+   stats.rst
    tensor/index.rst
+   rewriting.rst
    utilities/index.rst
-   parsing.rst
-   calculus/index.rst
-   physics/index.rst
-   categories.rst
-   diffgeom.rst
    vector/index.rst
-   algebras.rst
 
-Contributions to docs
+Contributions to Docs
 ---------------------
 
 All contributions are welcome. If you'd like to improve something, look into

doc/src/modules/integrals/integrals.rst

@@ -1,5 +1,6 @@
-Symbolic Integrals
-==================
+=========
+Integrals
+=========
 
 .. module:: sympy.integrals
 

doc/src/modules/liealgebras/index.rst

@@ -1,6 +1,6 @@
-==================
-Lie Algebra Module
-==================
+===========
+Lie Algebra
+===========
 
 .. automodule:: sympy.liealgebras
 

doc/src/modules/logic.rst

@@ -1,5 +1,6 @@
-Logic Module
-============
+=====
+Logic
+=====
 
 .. module:: sympy.logic
 

doc/src/modules/matrices/expressions.rst

@@ -34,6 +34,10 @@ Matrix Expressions Core Reference
    :members:
 .. autoclass:: MatPow
    :members:
+.. autoclass:: HadamardProduct
+   :members:
+.. autoclass:: HadamardPower
+   :members:
 .. autoclass:: Inverse
    :members:
 .. autoclass:: Transpose

doc/src/modules/ntheory.rst

@@ -1,5 +1,6 @@
+=============
 Number Theory
-====================
+=============
 
 .. module:: sympy.ntheory.generate
 
@@ -49,6 +50,8 @@ Ntheory Functions Reference
 
 .. autofunction:: factorint
 
+.. autofunction:: factorrat
+
 .. autofunction:: primefactors
 
 .. autofunction:: divisors
@@ -85,6 +88,20 @@ Ntheory Functions Reference
 .. autoclass:: primeomega
     :members:
 
+.. autofunction:: mersenne_prime_exponent
+
+.. autofunction:: is_perfect
+
+.. autofunction:: is_mersenne_prime
+
+.. autofunction:: abundance
+
+.. autofunction:: is_abundant
+
+.. autofunction:: is_deficient
+
+.. autofunction:: is_amicable
+
 .. module:: sympy.ntheory.modular
 
 .. autofunction:: symmetric_residue
@@ -113,8 +130,18 @@ Ntheory Functions Reference
 
 .. module:: sympy.ntheory.primetest
 
+.. autofunction:: is_euler_pseudoprime
+
+.. autofunction:: is_square
+
 .. autofunction:: mr
 
+.. autofunction:: is_lucas_prp
+
+.. autofunction:: is_strong_lucas_prp
+
+.. autofunction:: is_extra_strong_lucas_prp
+
 .. autofunction:: isprime
 
 .. module:: sympy.ntheory.residue_ntheory
@@ -127,6 +154,8 @@ Ntheory Functions Reference
 
 .. autofunction:: sqrt_mod
 
+.. autofunction:: sqrt_mod_iter
+
 .. autofunction:: quadratic_residues
 
 .. autofunction:: nthroot_mod
@@ -150,3 +179,7 @@ Ntheory Functions Reference
 .. module:: sympy.ntheory.egyptian_fraction
 
 .. autofunction:: egyptian_fraction
+
+.. module:: sympy.ntheory.bbp_pi
+
+.. autofunction:: pi_hex_digits

doc/src/modules/parsing.rst

@@ -1,5 +1,6 @@
-Parsing input
-=============
+=======
+Parsing
+=======
 
 Parsing Functions Reference
 ---------------------------
@@ -90,6 +91,27 @@ change between releases, and become stricter, more relaxed, or some mix.
 
 .. autoclass:: sympy.parsing.latex.LaTeXParsingError
 
+
+SymPy Expression Reference
+--------------------------
+
+.. module:: sympy.parsing.sym_expr
+
+.. autoclass:: SymPyExpression
+  :members:
+
+`Fortran` Parsing Reference
+---------------------------------
+
+.. module:: sympy.parsing.fortran.fortran_parser
+
+.. autoclass:: ASR2PyVisitor
+  :members:
+
+.. autofunction:: call_visitor
+
+.. autofunction:: src_to_sympy
+
 Runtime Installation
 --------------------
 

doc/src/modules/physics/continuum_mechanics/beam_problems.rst

@@ -175,7 +175,7 @@ final loading can be displayed:
        - 258⋅<x>   + 52⋅<x>   - 8⋅<x>  + 50⋅<x - 5>   + 8⋅<x - 5>  - 12⋅<x - 9>
 
 At this point, the beam is fully defined and the internal shear and bending
-moments are calculated like so:
+moments are calculated:
 
 .. plot::
    :context:
@@ -246,7 +246,7 @@ Example 2
 ---------
 
 There is a beam of length 30 meters. A moment of magnitude 120 Nm is
-applied in the clockwise direction at the end of the beam. A pointload
+applied in the clockwise direction at the end of the beam. A point load
 of magnitude 8 N is applied from the top of the beam at the starting
 point. There are two simple supports below the beam. One at the end
 and another one at a distance of 10 meters from the start. The
@@ -377,7 +377,7 @@ Example 4
 An overhanging beam of length 8 meters is pinned at 1 meter from starting point
 and supported by a roller 1 meter before the other end. It is subjected
 to a distributed constant load of 10 KN/m from the starting point till
-2 meters away from it. Two pointloads of 20KN and 8KN are applied at
+2 meters away from it. Two point loads of 20KN and 8KN are applied at
 5 meters and 7.5 meters away from the starting point respectively.
 
 ::
@@ -738,7 +738,7 @@ Example 9
 ---------
 
 There is a cantilever beam of length 4 meters. For first 2 meters
-its moment of inertia is ``1.5*I`` and ``I`` for the other end.
+its moment of inertia is ``1.5*I`` and ``I`` for the rest.
 A pointload of magnitude 20 N is applied from the top at its free end.
 
 ::
@@ -791,7 +791,7 @@ Example 10
 
 A combined beam, with constant flexural rigidity ``E*I``, is formed by joining
 a Beam of length ``2*l`` to the right of another Beam of length ``l``. The whole beam
-is fixed at both of its both end. A point load of magnitude ``P`` is also applied
+is fixed at both of its ends. A point load of magnitude ``P`` is also applied
 from the top at a distance of ``2*l`` from starting point.
 
 ::

doc/src/modules/physics/index.rst

@@ -1,8 +1,8 @@
 .. _physics-docs:
 
-==============
-Physics Module
-==============
+=======
+Physics
+=======
 
 .. automodule:: sympy.physics
 

doc/src/modules/physics/mechanics/autolev_parser.rst

@@ -82,7 +82,7 @@ The parser can be used as follows::
     >>> sympy_code = parse_autolev(open('double_pendulum.al'), include_numeric=True)
 
     # The include_pydy flag is False by default. Setting it to True will
-    # enable PyDy simulation code to be outputed if applicable.
+    # enable PyDy simulation code to be outputted if applicable.
 
     >>> print(sympy_code)
     import sympy.physics.mechanics as me

doc/src/modules/physics/mechanics/masses.rst

@@ -223,40 +223,42 @@ The following example shows how to use the momenta functions in
 One begins by creating the requisite symbols to describe the system. Then
 the reference frame is created and the kinematics are done. ::
 
-  >> from sympy import symbols
-  >> from sympy.physics.mechanics import dynamicsymbols, ReferenceFrame
-  >> from sympy.physics.mechanics import RigidBody, Particle, Point, outer
-  >> from sympy.physics.mechanics import linear_momentum, angular_momentum
-  >> m, M, l1 = symbols('m M l1')
-  >> q1d = dynamicsymbols('q1d')
-  >> N = ReferenceFrame('N')
-  >> O = Point('O')
-  >> O.set_vel(N, 0 * N.x)
-  >> Ac = O.locatenew('Ac', l1 * N.x)
-  >> P = Ac.locatenew('P', l1 * N.x)
-  >> a = ReferenceFrame('a')
-  >> a.set_ang_vel(N, q1d * N.z)
-  >> Ac.v2pt_theory(O, N, a)
-  >> P.v2pt_theory(O, N, a)
+  >>> from sympy import symbols
+  >>> from sympy.physics.mechanics import dynamicsymbols, ReferenceFrame
+  >>> from sympy.physics.mechanics import RigidBody, Particle, Point, outer
+  >>> from sympy.physics.mechanics import linear_momentum, angular_momentum
+  >>> m, M, l1 = symbols('m M l1')
+  >>> q1d = dynamicsymbols('q1d')
+  >>> N = ReferenceFrame('N')
+  >>> O = Point('O')
+  >>> O.set_vel(N, 0 * N.x)
+  >>> Ac = O.locatenew('Ac', l1 * N.x)
+  >>> P = Ac.locatenew('P', l1 * N.x)
+  >>> a = ReferenceFrame('a')
+  >>> a.set_ang_vel(N, q1d * N.z)
+  >>> Ac.v2pt_theory(O, N, a)
+  l1*q1d*N.y
+  >>> P.v2pt_theory(O, N, a)
+  2*l1*q1d*N.y
 
 Finally, the bodies that make up the system are created. In this case the
 system consists of a particle Pa and a RigidBody A. ::
 
-  >> Pa = Particle('Pa', P, m)
-  >> I = outer(N.z, N.z)
-  >> A = RigidBody('A', Ac, a, M, (I, Ac))
+  >>> Pa = Particle('Pa', P, m)
+  >>> I = outer(N.z, N.z)
+  >>> A = RigidBody('A', Ac, a, M, (I, Ac))
 
 Then one can either choose to evaluate the momenta of individual components
 of the system or of the entire system itself. ::
 
-  >> linear_momentum(N,A)
+  >>> linear_momentum(N,A)
   M*l1*q1d*N.y
-  >> angular_momentum(O, N, Pa)
+  >>> angular_momentum(O, N, Pa)
   4*l1**2*m*q1d*N.z
-  >> linear_momentum(N, A, Pa)
+  >>> linear_momentum(N, A, Pa)
   (M*l1*q1d + 2*l1*m*q1d)*N.y
-  >> angular_momentum(O, N, A, Pa)
-  (4*l1**2*m*q1d + q1d)*N.z
+  >>> angular_momentum(O, N, A, Pa)
+  (M*l1**2*q1d + 4*l1**2*m*q1d + q1d)*N.z
 
 It should be noted that the user can determine either momenta in any frame
 in :mod:`mechanics` as the user is allowed to specify the reference frame when
@@ -325,33 +327,34 @@ The following example shows how to use the energy functions in
 As was discussed above in the momenta functions, one first creates the system
 by going through an identical procedure. ::
 
-  >> from sympy import symbols
-  >> from sympy.physics.mechanics import dynamicsymbols, ReferenceFrame, outer
-  >> from sympy.physics.mechanics import RigidBody, Particle, mechanics_printing
-  >> from symp.physics.mechanics import kinetic_energy, potential_energy, Point
-  >> mechanics_printing()
-  >> m, M, l1, g, h, H = symbols('m M l1 g h H')
-  >> omega = dynamicsymbols('omega')
-  >> N = ReferenceFrame('N')
-  >> O = Point('O')
-  >> O.set_vel(N, 0 * N.x)
-  >> Ac = O.locatenew('Ac', l1 * N.x)
-  >> P = Ac.locatenew('P', l1 * N.x)
-  >> a = ReferenceFrame('a')
-  >> a.set_ang_vel(N, omega * N.z)
-  >> Ac.v2pt_theory(O, N, a)
-  >> P.v2pt_theory(O, N, a)
-  >> Pa = Particle('Pa', P, m)
-  >> I = outer(N.z, N.z)
-  >> A = RigidBody('A', Ac, a, M, (I, Ac))
+  >>> from sympy import symbols
+  >>> from sympy.physics.mechanics import dynamicsymbols, ReferenceFrame, outer
+  >>> from sympy.physics.mechanics import RigidBody, Particle
+  >>> from sympy.physics.mechanics import kinetic_energy, potential_energy, Point
+  >>> m, M, l1, g, h, H = symbols('m M l1 g h H')
+  >>> omega = dynamicsymbols('omega')
+  >>> N = ReferenceFrame('N')
+  >>> O = Point('O')
+  >>> O.set_vel(N, 0 * N.x)
+  >>> Ac = O.locatenew('Ac', l1 * N.x)
+  >>> P = Ac.locatenew('P', l1 * N.x)
+  >>> a = ReferenceFrame('a')
+  >>> a.set_ang_vel(N, omega * N.z)
+  >>> Ac.v2pt_theory(O, N, a)
+  l1*omega*N.y
+  >>> P.v2pt_theory(O, N, a)
+  2*l1*omega*N.y
+  >>> Pa = Particle('Pa', P, m)
+  >>> I = outer(N.z, N.z)
+  >>> A = RigidBody('A', Ac, a, M, (I, Ac))
 
 The user can then determine the kinetic energy of any number of entities of the
 system: ::
 
-  >> kinetic_energy(N, Pa)
-  2*l1**2*m*q1d**2
-  >> kinetic_energy(N, Pa, A)
-  M*l1**2*q1d**2/2 + 2*l1**2*m*q1d**2 + q1d**2/2
+  >>> kinetic_energy(N, Pa)
+  2*l1**2*m*omega(t)**2
+  >>> kinetic_energy(N, Pa, A)
+  M*l1**2*omega(t)**2/2 + 2*l1**2*m*omega(t)**2 + omega(t)**2/2
 
 It should be noted that the user can determine either kinetic energy relative
 to any frame in :mod:`mechanics` as the user is allowed to specify the
@@ -363,14 +366,15 @@ every entity of the system using the :mod:`potential_energy` property. The
 potential energy of any number of entities comprising the system can then be
 determined: ::
 
-  >> Pa.potential_energy = m * g * h
-  >> A.potential_energy = M * g * H
-  >> potential_energy(A, Pa)
+  >>> Pa.potential_energy = m * g * h
+  >>> A.potential_energy = M * g * H
+  >>> potential_energy(A, Pa)
   H*M*g + g*h*m
 
 One can also determine the Lagrangian for this system: ::
 
-  >> Lagrangian(Pa, A)
-  -H*M*g + M*l1**2*q1d**2/2 - g*h*m + 2*l1**2*m*q1d**2 + q1d**2/2
+  >>> from sympy.physics.mechanics import Lagrangian
+  >>> Lagrangian(N, Pa, A)
+  -H*M*g + M*l1**2*omega(t)**2/2 - g*h*m + 2*l1**2*m*omega(t)**2 + omega(t)**2/2
 
 Please refer to the docstrings to learn more about each function.

doc/src/modules/physics/units/examples.rst

@@ -108,9 +108,9 @@ the outcoming result):
 
     >>> from sympy.physics.units import convert_to
     >>> convert_to(q, day)
-    71.5123904642338*pi*day
+    71.5112118495813*pi*day
     >>> convert_to(q, day).n()
-    224.662800523082*day
+    224.659097795948*day
 
 We could also have the solar mass and the day as units coming from the
 astrophysical system, but we wanted to show how to create a unit that one needs.

doc/src/modules/plotting.rst

@@ -1,5 +1,6 @@
-Plotting Module
-===============
+========
+Plotting
+========
 
 .. module:: sympy.plotting.plot
 
@@ -87,8 +88,8 @@ Series Classes
 
 
 
-Pyglet Plotting Module
-======================
+Pyglet Plotting
+---------------
 
 .. module:: sympy.plotting.pygletplot
 

doc/src/modules/polys/index.rst

@@ -1,8 +1,8 @@
 .. _polys-docs:
 
-===============================
-Polynomials Manipulation Module
-===============================
+=======================
+Polynomial Manipulation
+=======================
 
 Computations with polynomials are at the core of computer algebra and
 having a fast and robust polynomials manipulation module is a key for

doc/src/modules/printing.rst

@@ -1,7 +1,8 @@
-Printing System
-===============
+========
+Printing
+========
 
-See the :ref:`tutorial-printing` section in Tutorial for introduction into
+See the :ref:`tutorial-printing` section in tutorial for introduction into
 printing.
 
 This guide documents the printing system in SymPy and how it works
@@ -13,8 +14,7 @@ Printer Class
 .. automodule:: sympy.printing.printer
 
 The main class responsible for printing is ``Printer`` (see also its
-`source code
-<https://github.com/sympy/sympy/blob/master/sympy/printing/printer.py>`_):
+`source code <https://github.com/sympy/sympy/blob/master/sympy/printing/printer.py>`_):
 
 .. autoclass:: Printer
     :members: doprint, _print, set_global_settings, order
@@ -185,10 +185,10 @@ An example where line wrapping is required:
 
     >>> expr = sqrt(1-x**2).series(x,n=20).removeO()
     >>> print(fcode(expr))
-          -715.0d0/65536.0d0*x**18 - 429.0d0/32768.0d0*x**16 - 33.0d0/
-         @ 2048.0d0*x**14 - 21.0d0/1024.0d0*x**12 - 7.0d0/256.0d0*x**10 -
-         @ 5.0d0/128.0d0*x**8 - 1.0d0/16.0d0*x**6 - 1.0d0/8.0d0*x**4 - 1.0d0
-         @ /2.0d0*x**2 + 1
+        -715.0d0/65536.0d0*x**18 - 429.0d0/32768.0d0*x**16 - 33.0d0/
+        @ 2048.0d0*x**14 - 21.0d0/1024.0d0*x**12 - 7.0d0/256.0d0*x**10 -
+        @ 5.0d0/128.0d0*x**8 - 1.0d0/16.0d0*x**6 - 1.0d0/8.0d0*x**4 - 1.0d0
+        @ /2.0d0*x**2 + 1
 
 In case of line wrapping, it is handy to include the assignment so that lines
 are wrapped properly when the assignment part is added.
@@ -301,6 +301,20 @@ Mathematica code printing
 
 .. autofunction:: sympy.printing.mathematica.mathematica_code
 
+Maple code printing
+-------------------
+
+.. module:: sympy.printing.maple
+
+.. autoclass:: sympy.printing.maple.MapleCodePrinter
+   :members:
+
+   .. autoattribute:: MapleCodePrinter.printmethod
+
+.. autofunction:: sympy.printing.maple.maple_code
+
+.. autofunction:: sympy.printing.maple.print_maple_code
+
 Javascript Code printing
 ------------------------
 

doc/src/modules/rewriting.rst

@@ -1,4 +1,5 @@
-Term rewriting
+==============
+Term Rewriting
 ==============
 
 Term rewriting is a very general class of functionalities which are used to
@@ -61,7 +62,7 @@ use ``rewrite()`` to do specific type of rewrite::
 Common Subexpression Detection and Collection
 ---------------------------------------------
 
-.. module:: sympy.simplify.cse_main
+.. currentmodule:: sympy.simplify.cse_main
 
 Before evaluating a large expression, it is often useful to identify common
 subexpressions, collect them and evaluate them at once. This is implemented

doc/src/modules/series/index.rst

@@ -1,8 +1,8 @@
 .. _series-docs:
 
-=============
-Series Module
-=============
+======
+Series
+======
 
 The series module implements series expansions as a function and many related
 functions.

doc/src/modules/sets.rst

@@ -1,5 +1,6 @@
+====
 Sets
-===========
+====
 
 .. automodule:: sympy.sets.sets
 
@@ -93,3 +94,27 @@ ComplexRegion
    :members:
 
 .. autofunction:: normalize_theta_set
+
+Iteration over sets
+^^^^^^^^^^^^^^^^^^^
+
+For set unions, `\{a, b\} \cup \{x, y\}` can be treated as
+`\{a, b, x, y\}` for iteration regardless of the distinctiveness of
+the elements, however, for set intersections, assuming that
+`\{a, b\} \cap \{x, y\}` is `\varnothing` or `\{a, b \}` would not
+always be valid, since some of `a`, `b`, `x` or `y` may or may not be
+the elements of the intersection.
+
+Iterating over the elements of a set involving intersection, complement,
+or symmetric difference yields (possibly duplicate) elements of the set
+provided that all elements are known to be the elements of the set.
+If any element cannot be determined to be a member of a set then the
+iteration gives ``TypeError``.
+This happens in the same cases where ``x in y`` would give an error.
+
+There are some reasons to implement like this, even if it breaks the
+consistency with how the python set iterator works.
+We keep in mind that sympy set comprehension like ``FiniteSet(*s)`` from
+a existing sympy sets could be a common usage.
+And this approach would make ``FiniteSet(*s)`` to be consistent with any
+symbolic set processing methods like ``FiniteSet(*simplify(s))``.

doc/src/modules/simplify/hyperexpand.rst

@@ -1,11 +1,12 @@
-Details on the Hypergeometric Function Expansion Module
-#######################################################
+========================
+Hypergeometric Expansion
+========================
 
 This page describes how the function :func:`hyperexpand` and related code
 work. For usage, see the documentation of the symplify module.
 
 Hypergeometric Function Expansion Algorithm
-*******************************************
+===========================================
 
 This section describes the algorithm used to expand hypergeometric functions.
 Most of it is based on the papers [Roach1996]_ and [Roach1997]_.

doc/src/modules/simplify/simplify.rst

@@ -1,5 +1,6 @@
+========
 Simplify
-********
+========
 
 .. currentmodule:: sympy.simplify.simplify
 
@@ -15,6 +16,10 @@ nthroot
 -------
 .. autofunction:: nthroot
 
+kroneckersimp
+-------------
+.. autofunction:: kroneckersimp
+
 besselsimp
 ----------
 .. autofunction:: besselsimp

doc/src/modules/solvers/ode.rst

@@ -64,6 +64,10 @@ the various ODE solving methods. For this reason, they are documented here.
 ^^^^^^^^^^^^^^^^^^^^^^^^^
 .. autofunction:: sympy.solvers.ode.ode_sol_simplicity
 
+:obj:`factorable`
+^^^^^^^^^^^^^^^^^^
+.. autofunction:: sympy.solvers.ode.ode_factorable
+
 :obj:`1st_exact`
 ^^^^^^^^^^^^^^^^
 .. autofunction:: sympy.solvers.ode.ode_1st_exact
@@ -84,6 +88,14 @@ the various ODE solving methods. For this reason, they are documented here.
 ^^^^^^^^^^^^^^^^^
 .. autofunction:: sympy.solvers.ode.ode_1st_linear
 
+:obj:`2nd_linear_airy`
+^^^^^^^^^^^^^^^^^^^^^^
+.. autofunction:: sympy.solvers.ode.ode_2nd_linear_airy
+
+:obj:`2nd_linear_bessel`
+^^^^^^^^^^^^^^^^^^^^^^^^
+.. autofunction:: sympy.solvers.ode.ode_2nd_linear_bessel
+
 :obj:`Bernoulli`
 ^^^^^^^^^^^^^^^^
 .. autofunction:: sympy.solvers.ode.ode_Bernoulli

doc/src/modules/solvers/solveset.rst

@@ -524,7 +524,7 @@ How are symbolic parameters handled in solveset?
     >>> not_empty_in(FiniteSet(x/2).intersect(Interval(0, 1)), x)
     Interval(0, 2)
     >>> not_empty_in(FiniteSet(x, x**2).intersect(Interval(1, 2)), x)
-    Union(Interval(-sqrt(2), -1), Interval(1, 2))
+    Union(Interval(1, 2), Interval(-sqrt(2), -1))
 
 
 References

doc/src/modules/stats.rst

@@ -1,5 +1,6 @@
+=====
 Stats
-===========
+=====
 
 .. automodule:: sympy.stats
 
@@ -7,7 +8,7 @@ Random Variable Types
 ^^^^^^^^^^^^^^^^^^^^^
 
 Finite Types
----------------
+------------
 .. autofunction:: DiscreteUniform
 .. autofunction:: Die
 .. autofunction:: Bernoulli
@@ -18,7 +19,7 @@ Finite Types
 .. autofunction:: Rademacher
 
 Discrete Types
------------------
+--------------
 .. autofunction:: Geometric
 .. autofunction:: Poisson
 .. autofunction:: Logarithmic
@@ -27,7 +28,7 @@ Discrete Types
 .. autofunction:: Zeta
 
 Continuous Types
--------------------
+----------------
 
 .. autofunction:: Arcsin
 .. autofunction:: Benini
@@ -69,6 +70,22 @@ Continuous Types
 .. autofunction:: WignerSemicircle
 .. autofunction:: ContinuousRV
 
+Joint Types
+-----------
+.. autofunction:: JointRV
+.. autofunction:: GeneralizedMultivariateLogGamma
+.. autofunction:: GeneralizedMultivariateLogGammaOmega
+.. autofunction:: Multinomial
+.. autofunction:: MultivariateBeta
+.. autofunction:: MultivariateEwens
+.. autofunction:: MultivariateT
+.. autofunction:: NegativeMultinomial
+.. autofunction:: NormalGamma
+
+Stochastic Processes
+--------------------
+.. autofunction:: DiscreteMarkovChain
+
 Interface
 ^^^^^^^^^
 
@@ -77,6 +94,7 @@ Interface
 .. autofunction:: E
 .. autoclass:: Expectation
 .. autofunction:: density
+.. autofunction:: entropy
 .. autofunction:: given
 .. autofunction:: where
 .. autofunction:: variance

doc/src/modules/tensor/index.rst

@@ -1,8 +1,8 @@
 .. _tensor_module:
 
-=============
-Tensor Module
-=============
+======
+Tensor
+======
 
 .. automodule:: sympy.tensor
 

doc/src/modules/tensor/tensor.rst

@@ -5,25 +5,12 @@ Tensor
 
 .. module:: sympy.tensor.tensor
 
-.. autoclass:: _TensorManager
-   :members:
-
 .. autoclass:: TensorIndexType
    :members:
 
 .. autoclass:: TensorIndex
    :members:
 
-.. autofunction:: tensor_indices
-
-.. autoclass:: TensorSymmetry
-   :members:
-
-.. autofunction:: tensorsymmetry
-
-.. autoclass:: TensorType
-   :members:
-
 .. autoclass:: TensorHead
    :members:
 
@@ -38,8 +25,17 @@ Tensor
 
 .. autofunction:: canon_bp
 
-.. autofunction:: tensor_mul
-
 .. autofunction:: riemann_cyclic_replace
 
 .. autofunction:: riemann_cyclic
+
+.. autoclass:: TensorSymmetry
+   :members:
+
+.. autofunction:: tensorsymmetry
+
+.. autoclass:: TensorType
+   :members:
+
+.. autoclass:: _TensorManager
+   :members:

doc/src/modules/vector/index.rst

@@ -1,6 +1,6 @@
-=============
-Vector Module
-=============
+======
+Vector
+======
 
 The vector module provides tools for basic vector math
 and differential calculus with respect to 3D Cartesian coordinate

doc/src/special_topics/finite_diff_derivatives.rst

@@ -101,7 +101,7 @@ Coefficients of `c_i` evaluated at `x_i + 2*h`:
 	>>> m32 = P(x0+2*h, x0, c, n).diff(c[1])
 	>>> m33 = P(x0+2*h, x0, c, n).diff(c[2])
 
-Matrix of the coeffcients is 3x3 in this case:
+Matrix of the coefficients is 3x3 in this case:
 
 	>>> M = Matrix([[m11, m12, m13], [m21, m22, m23], [m31, m32, m33]])
 
@@ -143,7 +143,7 @@ at `x_i` and thus using points at `x_{i-1}`,  `x_{i}`,  and `x_{i+1}`. So here i
 	>>> m31 = P(x0+h, x0, c, n).diff(c[0])
 	>>> m32 = P(x0+h, x0, c, n).diff(c[1])
 	>>> m33 = P(x0+h, x0, c, n).diff(c[2])
-	>>> # matrix of the coeffcients is 3x3 in this case
+	>>> # matrix of the coefficients is 3x3 in this case
 	>>> M = Matrix([[m11, m12, m13], [m21, m22, m23], [m31, m32, m33]])
 
 Now that we have the matrix of coefficients we next form the right-hand-side and solve by inverting `M`:
@@ -219,7 +219,7 @@ degree polynomial P coefficients of `c_i` evaluated at `x_i, x_{i-1},` and `x_{i
     >>> m32 = P(xN-2*h, xN, c, n).diff(c[1])
     >>> m33 = P(xN-2*h, xN, c, n).diff(c[2])
 
-Next we construct the `3 \times 3` matrix of the coeffcients:
+Next we construct the `3 \times 3` matrix of the coefficients:
 
     >>> M = Matrix([[m11, m12, m13], [m21, m22, m23], [m31, m32, m33]])
     >>> # matrix of the function values...actually a vector of right hand sides
@@ -247,7 +247,8 @@ Also,  we note that output of formats appropriate to Fortran,  C,  etc. may be d
 given above.
 
 Next we show how to perform these and many other discritizations of derivatives,  but using a
-much more efficient approach originally due to Bengt Fornberg and now incorported into SymPy.
-
+much more efficient approach originally due to Bengt Fornberg and now incorporated into SymPy.
 
+:ref:`calculus-finite-differences`
 
+:ref:`Finite difference weights <finite_diff>`

doc/src/tutorial/calculus.rst

@@ -196,11 +196,11 @@ definite integrals.  Here is a sampling of some of the power of ``integrate``.
     ⎮ sin⎝x ⎠ dx
     ⌡
     >>> integ.doit()
-                    ⎛√2⋅x⎞
-    3⋅√2⋅√π⋅fresnels⎜────⎟⋅Γ(3/4)
-                    ⎝ √π ⎠
-    ─────────────────────────────
-               8⋅Γ(7/4)
+             ⎛√2⋅x⎞
+    3⋅√2⋅√π⋅S⎜────⎟⋅Γ(3/4)
+             ⎝ √π ⎠
+    ──────────────────────
+           8⋅Γ(7/4)
 
     >>> integ = Integral(x**y*exp(-x), (x, 0, oo))
     >>> integ
@@ -322,6 +322,8 @@ The ``O`` notation supports arbitrary limit points (other than 0):
     -5 + ──────── + ──────── + ──────── + ──────── + x + O⎝(x - 6) ; x → 6⎠
             2          6          24        120
 
+.. _calculus-finite-differences:
+
 Finite differences
 ==================
 

doc/src/tutorial/manipulation.rst

@@ -30,37 +30,38 @@ expression looks like internally by using ``srepr``
 The easiest way to tear this apart is to look at a diagram of the expression
 tree:
 
-.. This comes from dotprint(2**x + x*y, labelfunc=srepr)
+.. This comes from dotprint(x**2 + x*y, labelfunc=srepr)
 
 .. graphviz::
 
     digraph{
 
     # Graph style
+    "ordering"="out"
     "rankdir"="TD"
 
     #########
     # Nodes #
     #########
 
-    "Symbol(x)_(0, 0)" ["color"="black", "label"="Symbol('x')", "shape"="ellipse"];
-    "Integer(2)_(1, 1)" ["color"="black", "label"="Integer(2)", "shape"="ellipse"];
-    "Symbol(y)_(0, 1)" ["color"="black", "label"="Symbol('y')", "shape"="ellipse"];
-    "Symbol(x)_(1, 0)" ["color"="black", "label"="Symbol('x')", "shape"="ellipse"];
-    "Mul(Symbol(x), Symbol(y))_(0,)" ["color"="black", "label"="Mul", "shape"="ellipse"];
-    "Pow(Symbol(x), Integer(2))_(1,)" ["color"="black", "label"="Pow", "shape"="ellipse"];
-    "Add(Mul(Symbol(x), Symbol(y)), Pow(Symbol(x), Integer(2)))_()" ["color"="black", "label"="Add", "shape"="ellipse"];
+    "Add(Pow(Symbol('x'), Integer(2)), Mul(Symbol('x'), Symbol('y')))_()" ["color"="black", "label"="Add", "shape"="ellipse"];
+    "Pow(Symbol('x'), Integer(2))_(0,)" ["color"="black", "label"="Pow", "shape"="ellipse"];
+    "Symbol('x')_(0, 0)" ["color"="black", "label"="Symbol('x')", "shape"="ellipse"];
+    "Integer(2)_(0, 1)" ["color"="black", "label"="Integer(2)", "shape"="ellipse"];
+    "Mul(Symbol('x'), Symbol('y'))_(1,)" ["color"="black", "label"="Mul", "shape"="ellipse"];
+    "Symbol('x')_(1, 0)" ["color"="black", "label"="Symbol('x')", "shape"="ellipse"];
+    "Symbol('y')_(1, 1)" ["color"="black", "label"="Symbol('y')", "shape"="ellipse"];
 
     #########
     # Edges #
     #########
 
-    "Mul(Symbol(x), Symbol(y))_(0,)" -> "Symbol(x)_(0, 0)";
-    "Mul(Symbol(x), Symbol(y))_(0,)" -> "Symbol(y)_(0, 1)";
-    "Pow(Symbol(x), Integer(2))_(1,)" -> "Symbol(x)_(1, 0)";
-    "Pow(Symbol(x), Integer(2))_(1,)" -> "Integer(2)_(1, 1)";
-    "Add(Mul(Symbol(x), Symbol(y)), Pow(Symbol(x), Integer(2)))_()" -> "Mul(Symbol(x), Symbol(y))_(0,)";
-    "Add(Mul(Symbol(x), Symbol(y)), Pow(Symbol(x), Integer(2)))_()" -> "Pow(Symbol(x), Integer(2))_(1,)";
+    "Add(Pow(Symbol('x'), Integer(2)), Mul(Symbol('x'), Symbol('y')))_()" -> "Pow(Symbol('x'), Integer(2))_(0,)";
+    "Add(Pow(Symbol('x'), Integer(2)), Mul(Symbol('x'), Symbol('y')))_()" -> "Mul(Symbol('x'), Symbol('y'))_(1,)";
+    "Pow(Symbol('x'), Integer(2))_(0,)" -> "Symbol('x')_(0, 0)";
+    "Pow(Symbol('x'), Integer(2))_(0,)" -> "Integer(2)_(0, 1)";
+    "Mul(Symbol('x'), Symbol('y'))_(1,)" -> "Symbol('x')_(1, 0)";
+    "Mul(Symbol('x'), Symbol('y'))_(1,)" -> "Symbol('y')_(1, 1)";
     }
 
 .. note::
@@ -82,16 +83,16 @@ expression, 2, we got ``Integer(2)``.  ``Integer`` is the SymPy class for
 integers.  It is similar to the Python built-in type ``int``, except that
 ``Integer`` plays nicely with other SymPy types.
 
-When we write ``2**x``, this creates a ``Pow`` object.  ``Pow`` is short for
+When we write ``x**2``, this creates a ``Pow`` object.  ``Pow`` is short for
 "power".
 
-    >>> srepr(2**x)
-    "Pow(Integer(2), Symbol('x'))"
+    >>> srepr(x**2)
+    "Pow(Symbol('x'), Integer(2))"
 
-We could have created the same object by calling ``Pow(2, x)``
+We could have created the same object by calling ``Pow(x, 2)``
 
-    >>> Pow(2, x)
-    2**x
+    >>> Pow(x, 2)
+    x**2
 
 Note that in the ``srepr`` output, we see ``Integer(2)``, the SymPy version of
 integers, even though technically, we input ``2``, a Python int.  In general,
@@ -104,7 +105,7 @@ function that does this is ``sympify`` [#sympify-fn]_.
     >>> type(sympify(2))
     <class 'sympy.core.numbers.Integer'>
 
-We have seen that ``2**x`` is represented as ``Pow(2, x)``.  What about
+We have seen that ``x**2`` is represented as ``Pow(x, 2)``.  What about
 ``x*y``?  As we might expect, this is the multiplication of ``x`` and ``y``.
 The SymPy class for multiplication is ``Mul``.
 
@@ -116,13 +117,13 @@ Thus, we could have created the same object by writing ``Mul(x, y)``.
     >>> Mul(x, y)
     x*y
 
-Now we get to our final expression, ``2**x + x*y``.  This is the addition of
-our last two objects, ``Pow(2, x)``, and ``Mul(x, y)``.  The SymPy class for
+Now we get to our final expression, ``x**2 + x*y``.  This is the addition of
+our last two objects, ``Pow(x, 2)``, and ``Mul(x, y)``.  The SymPy class for
 addition is ``Add``, so, as you might expect, to create this object, we use
-``Add(Pow(2, x), Mul(x, y))``.
+``Add(Pow(x, 2), Mul(x, y))``.
 
-    >>> Add(Pow(2, x), Mul(x, y))
-    2**x + x*y
+    >>> Add(Pow(x, 2), Mul(x, y))
+    x**2 + x*y
 
 SymPy expression trees can have many branches, and can be quite deep or quite
 broad.  Here is a more complicated example

doc/src/tutorial/matrices.rst

@@ -151,7 +151,7 @@ raise it to the ``-1`` power.
     >>> N**-1
     Traceback (most recent call last):
     ...
-    ValueError: Matrix det == 0; not invertible.
+    NonInvertibleMatrixError: Matrix det == 0; not invertible.
 
 To take the transpose of a Matrix, use ``T``.
 
@@ -427,7 +427,7 @@ guaranteed to be accurate in some limited domain of numerics and symbols,
 and any complicated expressions beyond its decidability are treated as ``None``,
 which behaves similarly to logical ``False``.
 
-The list of methods using zero testing procedures are as followings.
+The list of methods using zero testing procedures are as follows:
 
 ``echelon_form`` , ``is_echelon`` , ``rank`` , ``rref`` , ``nullspace`` ,
 ``eigenvects`` , ``inverse_ADJ`` , ``inverse_GE`` , ``inverse_LU`` ,

doc/src/tutorial/printing.rst

@@ -221,6 +221,7 @@ imported from ``sympy.printing.mathml``.
 ``print_mathml()`` prints the output.  If you want the string, use the
 function ``mathml()``.
 
+
 Dot
 ---
 
@@ -229,27 +230,33 @@ format, which can be rendered with Graphviz.  See the
 :ref:`tutorial-manipulation` section for some examples of the output of this
 printer.
 
+Here is an example of the raw output of the ``dotprint()`` function
+
     >>> from sympy.printing.dot import dotprint
     >>> from sympy.abc import x
-    >>> print(dotprint(x+2)) 
+    >>> print(dotprint(x+2))
     digraph{
+    <BLANKLINE>
     # Graph style
     "ordering"="out"
     "rankdir"="TD"
+    <BLANKLINE>
     #########
     # Nodes #
     #########
-    "Add(Integer(2), Symbol(x))_()" ["color"="black", "label"="Add", "shape"="ellipse"];
+    <BLANKLINE>
+    "Add(Integer(2), Symbol('x'))_()" ["color"="black", "label"="Add", "shape"="ellipse"];
     "Integer(2)_(0,)" ["color"="black", "label"="2", "shape"="ellipse"];
-    "Symbol(x)_(1,)" ["color"="black", "label"="x", "shape"="ellipse"];
+    "Symbol('x')_(1,)" ["color"="black", "label"="x", "shape"="ellipse"];
+    <BLANKLINE>
     #########
     # Edges #
     #########
-    "Add(Integer(2), Symbol(x))_()" -> "Integer(2)_(0,)";
-    "Add(Integer(2), Symbol(x))_()" -> "Symbol(x)_(1,)";
+    <BLANKLINE>
+    "Add(Integer(2), Symbol('x'))_()" -> "Integer(2)_(0,)";
+    "Add(Integer(2), Symbol('x'))_()" -> "Symbol('x')_(1,)";
     }
 
-
 .. rubric:: Footnotes
 
 .. [#srepr-fn] SymPy does not use the Python builtin ``repr()`` function for

doc/src/tutorial/solvers.rst

@@ -110,7 +110,7 @@ In the ``solveset`` module, the non linear system of equations is solved using
 	>>> nonlinsolve([x**2 + 1, y**2 + 1], [x, y])
 	{(-ⅈ, -ⅈ), (-ⅈ, ⅈ), (ⅈ, -ⅈ), (ⅈ, ⅈ)}
 
-3. When both real and complex solution is present:
+3. When both real and complex solution are present:
 
 	>>> from sympy import sqrt
 	>>> system = [x**2 - 2*y**2 -2, x*y - 2]
@@ -118,17 +118,11 @@ In the ``solveset`` module, the non linear system of equations is solved using
 	>>> nonlinsolve(system, vars)
 	{(-2, -1), (2, 1), (-√2⋅ⅈ, √2⋅ⅈ), (√2⋅ⅈ, -√2⋅ⅈ)}
 
-	>>> n = Dummy('n')
 	>>> system = [exp(x) - sin(y), 1/y - 3]
-	>>> real_soln = (log(sin(S(1)/3)), S(1)/3)
-	>>> img_lamda = Lambda(n, 2*n*I*pi + Mod(log(sin(S(1)/3)), 2*I*pi))
-	>>> complex_soln = (ImageSet(img_lamda, S.Integers), S(1)/3)
-	>>> soln = FiniteSet(real_soln, complex_soln)
-	>>> nonlinsolve(system, [x, y]) == soln
-	True
+	>>> nonlinsolve(system, vars)
+	{({2⋅n⋅ⅈ⋅π + log(sin(1/3)) | n ∊ ℤ}, 1/3)}
 
-4. If non linear system of equations is Positive dimensional system (A system with
-infinitely many solutions is said to be positive-dimensional):
+4. When the system is positive-dimensional system (has infinitely many solutions):
 
 	>>> nonlinsolve([x*y, x*y - x], [x, y])
 	{(0, y)}
@@ -148,14 +142,14 @@ infinitely many solutions is said to be positive-dimensional):
    ``solve`` can be used for such cases:
 
    >>> solve([x**2 - y**2/exp(x)], [x, y], dict=True)
-   ⎡⎧             ⎛y⎞⎫⎤
-   ⎢⎨x: 2⋅LambertW⎜─⎟⎬⎥
-   ⎣⎩             ⎝2⎠⎭⎦
+   ⎡⎧      ⎛-y ⎞⎫  ⎧      ⎛y⎞⎫⎤
+   ⎢⎨x: 2⋅W⎜───⎟⎬, ⎨x: 2⋅W⎜─⎟⎬⎥
+   ⎣⎩      ⎝ 2 ⎠⎭  ⎩      ⎝2⎠⎭⎦
 
    3. Currently ``nonlinsolve`` is not properly capable of solving the system of equations
    having trigonometric functions.
 
-   ``solve`` can be used for such cases(not all solution):
+   ``solve`` can be used for such cases (but does not give all solution):
 
    >>> solve([sin(x + y), cos(x - y)], [x, y])
    ⎡⎛-3⋅π   3⋅π⎞  ⎛-π   π⎞  ⎛π  3⋅π⎞  ⎛3⋅π  π⎞⎤
@@ -185,7 +179,7 @@ multiplicity 1 and ``3`` is a root of multiplicity 2.
    ``solve`` can be used for such cases:
 
    >>> solve(x*exp(x) - 1, x )
-   [LambertW(1)]
+   [W(1)]
 
 
 .. _tutorial-dsolve:

examples/advanced/pidigits.py

@@ -9,7 +9,7 @@ computation of the digits of pi.
 from mpmath import libmp, pi
 
 import math
-from time import clock
+from sympy.core.compatibility import clock
 import sys
 
 

examples/advanced/pyglet_plotting.py

@@ -8,10 +8,10 @@ Suggested Usage:    python -i pyglet_plotting.py
 
 
 from sympy import symbols, sin, cos, pi, sqrt
-from sympy.core.compatibility import range
+from sympy.core.compatibility import range, clock
 from sympy.plotting.pygletplot import PygletPlot
 
-from time import sleep, clock
+from time import sleep
 
 
 def main():

examples/beginner/limits_examples.py

@@ -23,7 +23,7 @@ def main():
 
     show( limit(x*(sqrt(x**2 + 1) - x), x, oo), Rational(1)/2 )
 
-    show( limit(x - sqrt(x**3 - 1), x, oo), Rational(0) )
+    show( limit(x - sqrt3(x**3 - 1), x, oo), Rational(0) )
 
     show( limit(log(1 + exp(x))/x, x, -oo), Rational(0) )
 

examples/intermediate/partial_differential_eqs.py

@@ -20,7 +20,7 @@ def main():
     eq = Eq(D(Xi(r, phi, theta), r, 2) + 2/r * D(Xi(r, phi, theta), r) +
             1/(r**2 * sin(phi)**2) * D(Xi(r, phi, theta), theta, 2) +
             cos(phi)/(r**2 * sin(phi)) * D(Xi(r, phi, theta), phi) +
-            1/r**2 * D(Xi(r, phi, theta), phi, 2))
+            1/r**2 * D(Xi(r, phi, theta), phi, 2), 0)
     pprint(eq)
 
     pprint("We can either separate this equation in regards with variable r:")

examples/intermediate/sample.py

@@ -21,7 +21,7 @@ def sample2d(f, x_args):
     try:
         f = sympify(f)
     except SympifyError:
-        raise ValueError("f could not be interpretted as a SymPy function")
+        raise ValueError("f could not be interpreted as a SymPy function")
     try:
         x, x_min, x_max, x_n = x_args
     except (TypeError, IndexError):

examples/notebooks/Macaulay_resultant.ipynb

@@ -44,7 +44,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "There are a number of steps needed to contruct matrix $A$. Let us consider an example from https://dl.acm.org/citation.cfm?id=550525 to \n",
+    "There are a number of steps needed to construct matrix $A$. Let us consider an example from https://dl.acm.org/citation.cfm?id=550525 to \n",
     "show the construction."
    ]
   },

examples/notebooks/Sylvester_resultant.ipynb

@@ -68,7 +68,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "From the definition, it is clear that the resultant will equal zero if and only if $p$ and $q$ have at least one common root. Thus, the resultant becomes very usefull in identifying whether common roots exist. "
+    "From the definition, it is clear that the resultant will equal zero if and only if $p$ and $q$ have at least one common root. Thus, the resultant becomes very useful in identifying whether common roots exist. "
    ]
   },
   {

release/rever.xsh

@@ -1110,8 +1110,6 @@ git_whitelist = {
     # This is just a distribute version of setup.py. Used mainly for setup.py
     # develop, which we don't care about in the release tarball
     'setupegg.py',
-    # Example on how to use tox to test Sympy. For development.
-    'tox.ini.sample',
     # pytest stuff
     'conftest.py',
     # Encrypted deploy key for deploying dev docs to GitHub

setup.py

@@ -33,8 +33,12 @@ import sys
 import os
 import shutil
 import glob
+import subprocess
 
-mpmath_version = '0.19'
+from distutils.command.sdist import sdist
+
+
+min_mpmath_version = '0.19'
 
 # This directory
 dir_setup = os.path.dirname(os.path.realpath(__file__))
@@ -58,11 +62,11 @@ except ImportError:
     from distutils.version import LooseVersion
     try:
         import mpmath
-        if mpmath.__version__ < LooseVersion(mpmath_version):
+        if mpmath.__version__ < LooseVersion(min_mpmath_version):
             raise ImportError
     except ImportError:
         print("Please install the mpmath package with a version >= %s"
-              % mpmath_version)
+              % min_mpmath_version)
         sys.exit(-1)
 
 PY3 = sys.version_info[0] > 2
@@ -122,6 +126,8 @@ modules = [
     'sympy.parsing.autolev._antlr',
     'sympy.parsing.autolev.test-examples',
     'sympy.parsing.autolev.test-examples.pydy-example-repo',
+    'sympy.parsing.c',
+    'sympy.parsing.fortran',
     'sympy.parsing.latex',
     'sympy.parsing.latex._antlr',
     'sympy.physics',
@@ -304,6 +310,37 @@ class antlr(Command):
             sys.exit(-1)
 
 
+class sdist_sympy(sdist):
+    def run(self):
+        # Fetch git commit hash and write down to commit_hash.txt before
+        # shipped in tarball.
+        commit_hash = None
+        commit_hash_filepath = 'doc/commit_hash.txt'
+        try:
+            commit_hash = \
+                subprocess.check_output(['git', 'rev-parse', 'HEAD'])
+            commit_hash = commit_hash.decode('ascii')
+            commit_hash = commit_hash.rstrip()
+            print('Commit hash found : {}.'.format(commit_hash))
+            print('Writing it to {}.'.format(commit_hash_filepath))
+        except:
+            pass
+
+        if commit_hash:
+            with open(commit_hash_filepath, 'w') as f:
+                f.write(commit_hash)
+
+        super(sdist_sympy, self).run()
+
+        try:
+            os.remove(commit_hash_filepath)
+            print(
+                'Successfully removed temporary file {}.'
+                .format(commit_hash_filepath))
+        except OSError as e:
+            print("Error deleting %s - %s." % (e.filename, e.strerror))
+
+
 # Check that this list is uptodate against the result of the command:
 # python bin/generate_test_list.py
 tests = [
@@ -409,6 +446,7 @@ if __name__ == '__main__':
                     'clean': clean,
                     'audit': audit,
                     'antlr': antlr,
+                    'sdist': sdist_sympy,
                     },
           classifiers=[
             'License :: OSI Approved :: BSD License',
@@ -420,14 +458,15 @@ if __name__ == '__main__':
             'Programming Language :: Python :: 2',
             'Programming Language :: Python :: 2.7',
             'Programming Language :: Python :: 3',
-            'Programming Language :: Python :: 3.4',
             'Programming Language :: Python :: 3.5',
             'Programming Language :: Python :: 3.6',
+            'Programming Language :: Python :: 3.7',
+            'Programming Language :: Python :: 3.8',
             'Programming Language :: Python :: Implementation :: CPython',
             'Programming Language :: Python :: Implementation :: PyPy',
             ],
           install_requires=[
-            'mpmath>=%s' % mpmath_version,
+            'mpmath>=%s' % min_mpmath_version,
             ],
           **extra_kwargs
           )

sympy/__init__.py

@@ -80,7 +80,7 @@ from .algebras import *
 # from combinatorics import *
 # This module is slow to import:
 #from physics import units
-from .plotting import plot, textplot, plot_backends, plot_implicit
+from .plotting import plot, textplot, plot_backends, plot_implicit, plot_parametric
 from .printing import *
 from .interactive import init_session, init_printing
 

sympy/abc.py

@@ -27,6 +27,28 @@ you still need to use ``Symbol('foo')`` or ``symbols('foo')``.
 You can freely mix usage of ``sympy.abc`` and ``Symbol``/``symbols``, though
 sticking with one and only one way to get the symbols does tend to make the code
 more readable.
+
+The module also defines some special names to help detect which names clash
+with the default SymPy namespace.
+
+``_clash1`` defines all the single letter variables that clash with
+SymPy objects; ``_clash2`` defines the multi-letter clashing symbols;
+and ``_clash`` is the union of both. These can be passed for ``locals``
+during sympification if one desires Symbols rather than the non-Symbol
+objects for those names.
+
+Examples
+========
+
+>>> from sympy import S
+>>> from sympy.abc import _clash1, _clash2, _clash
+>>> S("Q & C", locals=_clash1)
+C & Q
+>>> S('pi(x)', locals=_clash2)
+pi(x)
+>>> S('pi(C, Q)', locals=_clash)
+pi(C, Q)
+
 """
 
 from __future__ import print_function, division
@@ -70,50 +92,20 @@ _greek = list(greeks) # make a copy, so we can mutate it
 _greek.remove("lambda")
 _greek.append("lamda")
 
-def clashing():
-    """Return the clashing-symbols dictionaries.
+ns = {}
+exec_('from sympy import *', ns)
+_clash1 = {}
+_clash2 = {}
+while ns:
+    _k, _ = ns.popitem()
+    if _k in _greek:
+        _clash2[_k] = Symbol(_k)
+        _greek.remove(_k)
+    elif _k in _latin:
+        _clash1[_k] = Symbol(_k)
+        _latin.remove(_k)
+_clash = {}
+_clash.update(_clash1)
+_clash.update(_clash2)
 
-    ``clash1`` defines all the single letter variables that clash with
-    SymPy objects; ``clash2`` defines the multi-letter clashing symbols;
-    and ``clash`` is the union of both. These can be passed for ``locals``
-    during sympification if one desires Symbols rather than the non-Symbol
-    objects for those names.
-
-    Examples
-    ========
-
-    >>> from sympy import S
-    >>> from sympy.abc import _clash1, _clash2, _clash
-    >>> S("Q & C", locals=_clash1)
-    And(C, Q)
-    >>> S('pi(x)', locals=_clash2)
-    pi(x)
-    >>> S('pi(C, Q)', locals=_clash)
-    pi(C, Q)
-
-    Note: if changes are made to the docstring examples they can only
-    be tested after removing "clashing" from the list of deleted items
-    at the bottom of this file which removes this function from the
-    namespace.
-    """
-
-    ns = {}
-    exec_('from sympy import *', ns)
-    clash1 = {}
-    clash2 = {}
-    while ns:
-        k, _ = ns.popitem()
-        if k in _greek:
-            clash2[k] = Symbol(k)
-            _greek.remove(k)
-        elif k in _latin:
-            clash1[k] = Symbol(k)
-            _latin.remove(k)
-    clash = {}
-    clash.update(clash1)
-    clash.update(clash2)
-    return clash1, clash2, clash
-
-_clash1, _clash2, _clash = clashing()
-
-del _latin, _greek, clashing, Symbol
+del _latin, _greek, Symbol, _k

sympy/algebras/quaternion.py

@@ -3,16 +3,14 @@
 # https://en.wikipedia.org/wiki/Quaternion
 from __future__ import print_function
 
-from sympy import Rational
+from sympy import S, Rational
 from sympy import re, im, conjugate
 from sympy import sqrt, sin, cos, acos, exp, ln
 from sympy import trigsimp
 from sympy import integrate
-from sympy import Matrix, Add, Mul
+from sympy import Matrix
 from sympy import sympify
-from sympy.core.compatibility import SYMPY_INTS
 from sympy.core.expr import Expr
-from sympy.core.numbers import Integer
 
 
 class Quaternion(Expr):
@@ -110,8 +108,8 @@ class Quaternion(Expr):
         (x, y, z) = vector
         norm = sqrt(x**2 + y**2 + z**2)
         (x, y, z) = (x / norm, y / norm, z / norm)
-        s = sin(angle * Rational(1, 2))
-        a = cos(angle * Rational(1, 2))
+        s = sin(angle * S.Half)
+        a = cos(angle * S.Half)
         b = x * s
         c = y * s
         d = z * s
@@ -196,6 +194,16 @@ class Quaternion(Expr):
     def __neg__(self):
         return Quaternion(-self._a, -self._b, -self._c, -self.d)
 
+    def __truediv__(self, other):
+        return self * sympify(other)**-1
+
+    __div__ = __truediv__
+
+    def __rtruediv__(self, other):
+        return sympify(other) * self**-1
+
+    __rdiv__ = __rtruediv__
+
     def _eval_Integral(self, *args):
         return self.integrate(*args)
 
@@ -246,14 +254,12 @@ class Quaternion(Expr):
 
         # If q2 is a number or a sympy expression instead of a quaternion
         if not isinstance(q2, Quaternion):
-            if q1.real_field:
-                if q2.is_complex:
+            if q1.real_field and q2.is_complex:
                     return Quaternion(re(q2) + q1.a, im(q2) + q1.b, q1.c, q1.d)
-                else:
-                    # q2 is something strange, do not evaluate:
-                    return Add(q1, q2)
-            else:
+            elif q2.is_commutative:
                 return Quaternion(q1.a + q2, q1.b, q1.c, q1.d)
+            else:
+                raise ValueError("Only commutative expressions can be added with a Quaternion.")
 
         return Quaternion(q1.a + q2.a, q1.b + q2.b, q1.c + q2.c, q1.d
                           + q2.d)
@@ -321,14 +327,14 @@ class Quaternion(Expr):
         ========
 
         >>> from sympy.algebras.quaternion import Quaternion
-        >>> from sympy import symbols
+        >>> from sympy import Symbol
         >>> q1 = Quaternion(1, 2, 3, 4)
         >>> q2 = Quaternion(5, 6, 7, 8)
         >>> Quaternion._generic_mul(q1, q2)
         (-60) + 12*i + 30*j + 24*k
         >>> Quaternion._generic_mul(q1, 2)
         2 + 4*i + 6*j + 8*k
-        >>> x = symbols('x', real = True)
+        >>> x = Symbol('x', real = True)
         >>> Quaternion._generic_mul(q1, x)
         x + 2*x*i + 3*x*j + 4*x*k
 
@@ -349,24 +355,21 @@ class Quaternion(Expr):
 
         # If q1 is a number or a sympy expression instead of a quaternion
         if not isinstance(q1, Quaternion):
-            if q2.real_field:
-                if q1.is_complex:
-                    return q2 * Quaternion(re(q1), im(q1), 0, 0)
-                else:
-                    return Mul(q1, q2)
-            else:
+            if q2.real_field and q1.is_complex:
+                    return Quaternion(re(q1), im(q1), 0, 0) * q2
+            elif q1.is_commutative:
                 return Quaternion(q1 * q2.a, q1 * q2.b, q1 * q2.c, q1 * q2.d)
-
+            else:
+                raise ValueError("Only commutative expressions can be multiplied with a Quaternion.")
 
         # If q2 is a number or a sympy expression instead of a quaternion
         if not isinstance(q2, Quaternion):
-            if q1.real_field:
-                if q2.is_complex:
+            if q1.real_field and q2.is_complex:
                     return q1 * Quaternion(re(q2), im(q2), 0, 0)
-                else:
-                    return Mul(q1, q2)
-            else:
+            elif q2.is_commutative:
                 return Quaternion(q2 * q1.a, q2 * q1.b, q2 * q1.c, q2 * q1.d)
+            else:
+                raise ValueError("Only commutative expressions can be multiplied with a Quaternion.")
 
         return Quaternion(-q1.b*q2.b - q1.c*q2.c - q1.d*q2.d + q1.a*q2.a,
                           q1.b*q2.a + q1.c*q2.d - q1.d*q2.c + q1.a*q2.b,
@@ -421,22 +424,23 @@ class Quaternion(Expr):
         >>> q.pow(4)
         668 + (-224)*i + (-336)*j + (-448)*k
         """
+        p = sympify(p)
         q = self
         if p == -1:
             return q.inverse()
         res = 1
 
+        if not p.is_Integer:
+            return NotImplemented
+
         if p < 0:
             q, p = q.inverse(), -p
 
-        if not (isinstance(p, (Integer, SYMPY_INTS))):
-            return NotImplemented
-
         while p > 0:
-            if p & 1:
+            if p % 2 == 1:
                 res = q * res
 
-            p = p >> 1
+            p = p//2
             q = q * q
 
         return res
@@ -597,14 +601,8 @@ class Quaternion(Expr):
         2*pi/3
         """
         q = self
-        try:
-            # Skips it if it doesn't know whether q.a is negative
-            if q.a < 0:
-                # avoid error with acos
-                # axis and angle of rotation of q and q*-1 will be the same
-                q = q * -1
-        except BaseException:
-            pass
+        if q.a.is_negative:
+            q = q * -1
 
         q = q.normalize()
         angle = trigsimp(2 * acos(q.a))

sympy/algebras/tests/test_quaternion.py

@@ -1,10 +1,8 @@
+from sympy import symbols, re, im, I, Abs, Symbol, \
+     cos, sin, sqrt, conjugate, log, acos, E, pi, \
+     Matrix, diff, integrate, trigsimp, S, Rational
 from sympy.algebras.quaternion import Quaternion
-from sympy import symbols, re, im, Add, Mul, I, Abs
-from sympy import cos, sin, sqrt, conjugate, log, acos, E, pi
 from sympy.utilities.pytest import raises
-from sympy import Matrix
-from sympy import diff, integrate, trigsimp
-from sympy import S, Rational
 
 x, y, z, w = symbols("x y z w")
 
@@ -12,14 +10,18 @@ def test_quaternion_construction():
     q = Quaternion(x, y, z, w)
     assert q + q == Quaternion(2*x, 2*y, 2*z, 2*w)
 
-    q2 = Quaternion.from_axis_angle((sqrt(3)/3, sqrt(3)/3, sqrt(3)/3), 2*pi/3)
-    assert q2 == Quaternion(Rational(1, 2), Rational(1, 2),
-                            Rational(1, 2), Rational(1, 2))
+    q2 = Quaternion.from_axis_angle((sqrt(3)/3, sqrt(3)/3, sqrt(3)/3),
+                                    pi*Rational(2, 3))
+    assert q2 == Quaternion(S.Half, S.Half,
+                            S.Half, S.Half)
 
     M = Matrix([[cos(x), -sin(x), 0], [sin(x), cos(x), 0], [0, 0, 1]])
     q3 = trigsimp(Quaternion.from_rotation_matrix(M))
     assert q3 == Quaternion(sqrt(2)*sqrt(cos(x) + 1)/2, 0, 0, sqrt(-2*cos(x) + 2)/2)
 
+    nc = Symbol('nc', commutative=False)
+    raises(ValueError, lambda: Quaternion(x, y, nc, w))
+
 
 def test_quaternion_complex_real_addition():
     a = symbols("a", complex=True)
@@ -32,9 +34,9 @@ def test_quaternion_complex_real_addition():
     assert 1 + q == Quaternion(1 + x, y, z, w)
     assert I + q == Quaternion(x, 1 + y, z, w)
     assert b + q == Quaternion(x + b, y, z, w)
-    assert c + q == Add(c, Quaternion(x, y, z, w), evaluate=False)
-    assert q * c == Mul(Quaternion(x, y, z, w), c, evaluate=False)
-    assert c * q == Mul(c, Quaternion(x, y, z, w), evaluate=False)
+    raises(ValueError, lambda: c + q)
+    raises(ValueError, lambda: q * c)
+    raises(ValueError, lambda: c * q)
 
     assert -q == Quaternion(-x, -y, -z, -w)
 
@@ -47,6 +49,12 @@ def test_quaternion_complex_real_addition():
     Quaternion((2 + 3*I)*(3 + 4*I), (2 + 3*I)*(2 + 5*I), 0, (2 + 3*I)*(7 + 8*I))
     assert q2 * (2 + 3*I) == Quaternion(-10, 11, 38, -5)
 
+    q1 = Quaternion(1, 2, 3, 4)
+    q0 = Quaternion(0, 0, 0, 0)
+    assert q1 + q0 == q1
+    assert q1 - q0 == q1
+    assert q1 - q1 == q0
+
 
 def test_quaternion_functions():
     q = Quaternion(x, y, z, w)
@@ -57,10 +65,14 @@ def test_quaternion_functions():
     assert q.norm() == sqrt(w**2 + x**2 + y**2 + z**2)
     assert q.normalize() == Quaternion(x, y, z, w) / sqrt(w**2 + x**2 + y**2 + z**2)
     assert q.inverse() == Quaternion(x, -y, -z, -w) / (w**2 + x**2 + y**2 + z**2)
+    assert q.inverse() == q.pow(-1)
     raises(ValueError, lambda: q0.inverse())
     assert q.pow(2) == Quaternion(-w**2 + x**2 - y**2 - z**2, 2*x*y, 2*x*z, 2*w*x)
-    assert q1.pow(-2) == Quaternion(-S(7)/225, -S(1)/225, -S(1)/150, -S(2)/225)
+    assert q**(2) == Quaternion(-w**2 + x**2 - y**2 - z**2, 2*x*y, 2*x*z, 2*w*x)
+    assert q1.pow(-2) == Quaternion(Rational(-7, 225), Rational(-1, 225), Rational(-1, 150), Rational(-2, 225))
+    assert q1**(-2) == Quaternion(Rational(-7, 225), Rational(-1, 225), Rational(-1, 150), Rational(-2, 225))
     assert q1.pow(-0.5) == NotImplemented
+    raises(TypeError, lambda: q1**(-0.5))
 
     assert q1.exp() == \
     Quaternion(E * cos(sqrt(29)),
@@ -84,7 +96,11 @@ def test_quaternion_functions():
     assert integrate(Quaternion(x, x, x, x), x) == \
     Quaternion(x**2 / 2, x**2 / 2, x**2 / 2, x**2 / 2)
 
-    assert Quaternion.rotate_point((1, 1, 1), q1) == (S(1) / 5, 1, S(7) / 5)
+    assert Quaternion.rotate_point((1, 1, 1), q1) == (S.One / 5, 1, S(7) / 5)
+    n = Symbol('n')
+    raises(TypeError, lambda: q1**n)
+    n = Symbol('n', integer=True)
+    raises(TypeError, lambda: q1**n)
 
 
 def test_quaternion_conversions():
@@ -95,14 +111,14 @@ def test_quaternion_conversions():
                                    4 * sqrt(29)/29),
                                    2 * acos(sqrt(30)/30))
 
-    assert q1.to_rotation_matrix() == Matrix([[-S(2)/3, S(2)/15, S(11)/15],
-                                     [S(2)/3, -S(1)/3, S(2)/3],
-                                     [S(1)/3, S(14)/15, S(2)/15]])
+    assert q1.to_rotation_matrix() == Matrix([[Rational(-2, 3), Rational(2, 15), Rational(11, 15)],
+                                              [Rational(2, 3), Rational(-1, 3), Rational(2, 3)],
+                                              [Rational(1, 3), Rational(14, 15), Rational(2, 15)]])
 
-    assert q1.to_rotation_matrix((1, 1, 1)) == Matrix([[-S(2)/3, S(2)/15, S(11)/15, S(4)/5],
-                                                  [S(2)/3, -S(1)/3, S(2)/3, S(0)],
-                                                       [S(1)/3, S(14)/15, S(2)/15, -S(2)/5],
-                                                  [S(0), S(0), S(0), S(1)]])
+    assert q1.to_rotation_matrix((1, 1, 1)) == Matrix([[Rational(-2, 3), Rational(2, 15), Rational(11, 15), Rational(4, 5)],
+                                                       [Rational(2, 3), Rational(-1, 3), Rational(2, 3), S.Zero],
+                                                       [Rational(1, 3), Rational(14, 15), Rational(2, 15), Rational(-2, 5)],
+                                                       [S.Zero, S.Zero, S.Zero, S.One]])
 
     theta = symbols("theta", real=True)
     q2 = Quaternion(cos(theta/2), 0, 0, sin(theta/2))
@@ -136,3 +152,22 @@ def test_quaternion_rotation_iss1593():
                 [1,      0,      0],
                 [0, cos(x), -sin(x)],
                 [0, sin(x), cos(x)]]))
+
+
+def test_quaternion_multiplication():
+    q1 = Quaternion(3 + 4*I, 2 + 5*I, 0, 7 + 8*I, real_field = False)
+    q2 = Quaternion(1, 2, 3, 5)
+    q3 = Quaternion(1, 1, 1, y)
+
+    assert Quaternion._generic_mul(4, 1) == 4
+    assert Quaternion._generic_mul(4, q1) == Quaternion(12 + 16*I, 8 + 20*I, 0, 28 + 32*I)
+    assert q2.mul(2) == Quaternion(2, 4, 6, 10)
+    assert q2.mul(q3) == Quaternion(-5*y - 4, 3*y - 2, 9 - 2*y, y + 4)
+    assert q2.mul(q3) == q2*q3
+
+    z = symbols('z', complex=True)
+    z_quat = Quaternion(re(z), im(z), 0, 0)
+    q = Quaternion(*symbols('q:4', real=True))
+
+    assert z * q == z_quat * q
+    assert q * z == q * z_quat

sympy/assumptions/ask.py

@@ -8,9 +8,10 @@ from sympy.core.cache import cacheit
 from sympy.core.decorators import deprecated
 from sympy.core.relational import Relational
 from sympy.logic.boolalg import (to_cnf, And, Not, Or, Implies, Equivalent,
-    BooleanFunction, BooleanAtom)
+                                 BooleanFunction, BooleanAtom)
 from sympy.logic.inference import satisfiable
 from sympy.utilities.decorator import memoize_property
+from sympy.assumptions.cnf import CNF, EncodedCNF, Literal
 
 
 # Deprecated predicates should be added to this list
@@ -528,7 +529,7 @@ class AssumptionKeys(object):
         False
         >>> ask(~Q.zero(I))
         True
-        >>> ask(Q.nonzero(oo))  #doctest: +SKIP
+        >>> ask(Q.nonzero(oo))
         False
 
         """
@@ -553,6 +554,7 @@ class AssumptionKeys(object):
         ========
 
         >>> from sympy import Q, ask, I
+
         >>> ask(Q.nonpositive(-1))
         True
         >>> ask(Q.nonpositive(0))
@@ -1211,6 +1213,28 @@ def _extract_facts(expr, symbol, check_reversed_rel=True):
         return expr.func(*args)
 
 
+def _extract_all_facts(expr, symbol):
+    facts = set()
+    if isinstance(symbol, Relational):
+        symbols = (symbol, symbol.reversed)
+    else:
+        symbols = (symbol,)
+    for clause in expr.clauses:
+        args = []
+        for literal in clause:
+            if isinstance(literal.lit, AppliedPredicate):
+                if literal.lit.arg in symbols:
+                    # Add literal if it has 'symbol' in it
+                    args.append(Literal(literal.lit.func, literal.is_Not))
+                else:
+                    # If any of the literals doesn't have 'symbol' don't add the whole clause.
+                    break
+        else:
+            if args:
+                facts.add(frozenset(args))
+    return CNF(facts)
+
+
 def ask(proposition, assumptions=True, context=global_assumptions):
     """
     Method for inferring properties about objects.
@@ -1239,7 +1263,7 @@ def ask(proposition, assumptions=True, context=global_assumptions):
         Relations in assumptions are not implemented (yet), so the following
         will not give a meaningful result.
 
-        >>> ask(Q.positive(x), Q.is_true(x > 0)) # doctest: +SKIP
+        >>> ask(Q.positive(x), Q.is_true(x > 0))
 
         It is however a work in progress.
 
@@ -1257,54 +1281,54 @@ def ask(proposition, assumptions=True, context=global_assumptions):
     else:
         key, expr = Q.is_true, sympify(proposition)
 
-    assumptions = And(assumptions, And(*context))
-    assumptions = to_cnf(assumptions)
+    assump = CNF.from_prop(assumptions)
+    assump.extend(context)
 
-    local_facts = _extract_facts(assumptions, expr)
+    local_facts = _extract_all_facts(assump, expr)
 
-    known_facts_cnf = get_known_facts_cnf()
+    known_facts_cnf = get_all_known_facts()
     known_facts_dict = get_known_facts_dict()
 
-    if local_facts and satisfiable(And(local_facts, known_facts_cnf)) is False:
+    enc_cnf = EncodedCNF()
+    enc_cnf.from_cnf(CNF(known_facts_cnf))
+    enc_cnf.add_from_cnf(local_facts)
+
+    if local_facts.clauses and satisfiable(enc_cnf) is False:
         raise ValueError("inconsistent assumptions %s" % assumptions)
 
+    if local_facts.clauses:
+        local_facts_ = CNF.CNF_to_cnf(local_facts)
+
+        # See if there's a straight-forward conclusion we can make for the inference
+        if local_facts_.is_Atom:
+            if key in known_facts_dict[local_facts_]:
+                return True
+            if Not(key) in known_facts_dict[local_facts_]:
+                return False
+        elif (isinstance(local_facts_, And) and
+              all(k in known_facts_dict for k in local_facts_.args)):
+            for assum in local_facts_.args:
+                if assum.is_Atom:
+                    if key in known_facts_dict[assum]:
+                        return True
+                    if Not(key) in known_facts_dict[assum]:
+                        return False
+                elif isinstance(assum, Not) and assum.args[0].is_Atom:
+                    if key in known_facts_dict[assum]:
+                        return False
+                    if Not(key) in known_facts_dict[assum]:
+                        return True
+        elif (isinstance(key, Predicate) and
+              isinstance(local_facts_, Not) and local_facts_.args[0].is_Atom):
+            if local_facts_.args[0] in known_facts_dict[key]:
+                return False
+
     # direct resolution method, no logic
     res = key(expr)._eval_ask(assumptions)
     if res is not None:
         return bool(res)
-
-    if local_facts is None:
-        return satask(proposition, assumptions=assumptions, context=context)
-
-
-    # See if there's a straight-forward conclusion we can make for the inference
-    if local_facts.is_Atom:
-        if key in known_facts_dict[local_facts]:
-            return True
-        if Not(key) in known_facts_dict[local_facts]:
-            return False
-    elif (isinstance(local_facts, And) and
-            all(k in known_facts_dict for k in local_facts.args)):
-        for assum in local_facts.args:
-            if assum.is_Atom:
-                if key in known_facts_dict[assum]:
-                    return True
-                if Not(key) in known_facts_dict[assum]:
-                    return False
-            elif isinstance(assum, Not) and assum.args[0].is_Atom:
-                if key in known_facts_dict[assum]:
-                    return False
-                if Not(key) in known_facts_dict[assum]:
-                    return True
-    elif (isinstance(key, Predicate) and
-            isinstance(local_facts, Not) and local_facts.args[0].is_Atom):
-        if local_facts.args[0] in known_facts_dict[key]:
-            return False
-
-    # Failing all else, we do a full logical inference
-    res = ask_full_inference(key, local_facts, known_facts_cnf)
-    if res is None:
-        return satask(proposition, assumptions=assumptions, context=context)
+    # using satask (still costly)
+    res = satask(proposition, assumptions=assumptions, context=context)
     return res
 
 
@@ -1387,8 +1411,16 @@ def compute_known_facts(known_facts, known_facts_keys):
 
     from sympy.core.cache import cacheit
     from sympy.logic.boolalg import And
+    from sympy.assumptions.cnf import Literal
     from sympy.assumptions.ask import Q
 
+    # -{ Known facts as a set }-
+    @cacheit
+    def get_all_known_facts():
+        return {
+            %s
+        }
+
     # -{ Known facts in Conjunctive Normal Form }-
     @cacheit
     def get_known_facts_cnf():
@@ -1407,7 +1439,10 @@ def compute_known_facts(known_facts, known_facts_keys):
     LINE = ",\n        "
     HANG = ' '*8
     cnf = to_cnf(known_facts)
+    cnf_ = CNF.to_CNF(known_facts)
     c = LINE.join([str(a) for a in cnf.args])
+
+    p = LINE.join(sorted(['frozenset((' + ', '.join(str(lit) for lit in sorted(clause, key=str)) +'))' for clause in cnf_.clauses]))
     mapping = single_fact_lookup(known_facts_keys, cnf)
     items = sorted(mapping.items(), key=str)
     keys = [str(i[0]) for i in items]
@@ -1417,7 +1452,7 @@ def compute_known_facts(known_facts, known_facts_keys):
             subsequent_indent=HANG,
             break_long_words=False))
         for k, v in zip(keys, values)]) + ','
-    return fact_string % (c, m)
+    return fact_string % (p, c, m)
 
 # handlers tells us what ask handler we should use
 # for a particular key
@@ -1481,10 +1516,9 @@ def get_known_facts():
         Equivalent(Q.extended_real, Q.real | Q.infinite),
         Equivalent(Q.even | Q.odd, Q.integer),
         Implies(Q.even, ~Q.odd),
-        Equivalent(Q.prime, Q.integer & Q.positive & ~Q.composite),
+        Implies(Q.prime, Q.integer & Q.positive & ~Q.composite),
         Implies(Q.integer, Q.rational),
         Implies(Q.rational, Q.algebraic),
-        Implies(Q.irrational, Q.finite),
         Implies(Q.algebraic, Q.complex),
         Implies(Q.algebraic, Q.finite),
         Equivalent(Q.transcendental | Q.algebraic, Q.complex & Q.finite),
@@ -1530,4 +1564,4 @@ def get_known_facts():
     )
 
 from sympy.assumptions.ask_generated import (
-    get_known_facts_dict, get_known_facts_cnf)
+    get_known_facts_dict, get_known_facts_cnf, get_all_known_facts)

sympy/assumptions/ask_generated.py

@@ -8,8 +8,87 @@ Instead, run ./bin/ask_update.py.
 
 from sympy.core.cache import cacheit
 from sympy.logic.boolalg import And
+from sympy.assumptions.cnf import Literal
 from sympy.assumptions.ask import Q
 
+# -{ Known facts as a set }-
+@cacheit
+def get_all_known_facts():
+    return {
+        frozenset((Literal(Q.algebraic, False), Literal(Q.complex, True), Literal(Q.finite, True), Literal(Q.transcendental, False))),
+        frozenset((Literal(Q.algebraic, False), Literal(Q.rational, True))),
+        frozenset((Literal(Q.algebraic, True), Literal(Q.complex, False))),
+        frozenset((Literal(Q.algebraic, True), Literal(Q.finite, False))),
+        frozenset((Literal(Q.algebraic, True), Literal(Q.transcendental, True))),
+        frozenset((Literal(Q.antihermitian, False), Literal(Q.imaginary, True))),
+        frozenset((Literal(Q.antihermitian, True), Literal(Q.hermitian, True))),
+        frozenset((Literal(Q.complex, False), Literal(Q.imaginary, True))),
+        frozenset((Literal(Q.complex, False), Literal(Q.real, True))),
+        frozenset((Literal(Q.complex, False), Literal(Q.transcendental, True))),
+        frozenset((Literal(Q.complex_elements, False), Literal(Q.real_elements, True))),
+        frozenset((Literal(Q.composite, True), Literal(Q.prime, True))),
+        frozenset((Literal(Q.diagonal, False), Literal(Q.lower_triangular, True), Literal(Q.upper_triangular, True))),
+        frozenset((Literal(Q.diagonal, True), Literal(Q.lower_triangular, False))),
+        frozenset((Literal(Q.diagonal, True), Literal(Q.normal, False))),
+        frozenset((Literal(Q.diagonal, True), Literal(Q.symmetric, False))),
+        frozenset((Literal(Q.diagonal, True), Literal(Q.upper_triangular, False))),
+        frozenset((Literal(Q.even, False), Literal(Q.integer, True), Literal(Q.odd, False))),
+        frozenset((Literal(Q.even, False), Literal(Q.zero, True))),
+        frozenset((Literal(Q.even, True), Literal(Q.integer, False))),
+        frozenset((Literal(Q.even, True), Literal(Q.odd, True))),
+        frozenset((Literal(Q.extended_real, False), Literal(Q.infinite, True))),
+        frozenset((Literal(Q.extended_real, False), Literal(Q.real, True))),
+        frozenset((Literal(Q.extended_real, True), Literal(Q.infinite, False), Literal(Q.real, False))),
+        frozenset((Literal(Q.finite, False), Literal(Q.irrational, True))),
+        frozenset((Literal(Q.finite, False), Literal(Q.rational, True))),
+        frozenset((Literal(Q.finite, False), Literal(Q.transcendental, True))),
+        frozenset((Literal(Q.finite, True), Literal(Q.infinite, True))),
+        frozenset((Literal(Q.finite, True), Literal(Q.irrational, False), Literal(Q.rational, False), Literal(Q.real, True))),
+        frozenset((Literal(Q.fullrank, False), Literal(Q.invertible, True))),
+        frozenset((Literal(Q.fullrank, True), Literal(Q.invertible, False), Literal(Q.square, True))),
+        frozenset((Literal(Q.hermitian, False), Literal(Q.real, True))),
+        frozenset((Literal(Q.imaginary, True), Literal(Q.real, True))),
+        frozenset((Literal(Q.integer, False), Literal(Q.odd, True))),
+        frozenset((Literal(Q.integer, False), Literal(Q.prime, True))),
+        frozenset((Literal(Q.integer, True), Literal(Q.rational, False))),
+        frozenset((Literal(Q.integer_elements, True), Literal(Q.real_elements, False))),
+        frozenset((Literal(Q.invertible, False), Literal(Q.positive_definite, True))),
+        frozenset((Literal(Q.invertible, False), Literal(Q.singular, False))),
+        frozenset((Literal(Q.invertible, False), Literal(Q.unitary, True))),
+        frozenset((Literal(Q.invertible, True), Literal(Q.singular, True))),
+        frozenset((Literal(Q.invertible, True), Literal(Q.square, False))),
+        frozenset((Literal(Q.irrational, True), Literal(Q.rational, True))),
+        frozenset((Literal(Q.irrational, True), Literal(Q.real, False))),
+        frozenset((Literal(Q.lower_triangular, False), Literal(Q.triangular, True), Literal(Q.upper_triangular, False))),
+        frozenset((Literal(Q.lower_triangular, True), Literal(Q.triangular, False))),
+        frozenset((Literal(Q.negative, False), Literal(Q.nonpositive, True), Literal(Q.zero, False))),
+        frozenset((Literal(Q.negative, False), Literal(Q.nonzero, True), Literal(Q.positive, False))),
+        frozenset((Literal(Q.negative, False), Literal(Q.positive, False), Literal(Q.real, True), Literal(Q.zero, False))),
+        frozenset((Literal(Q.negative, True), Literal(Q.nonpositive, False))),
+        frozenset((Literal(Q.negative, True), Literal(Q.nonzero, False))),
+        frozenset((Literal(Q.negative, True), Literal(Q.positive, True))),
+        frozenset((Literal(Q.negative, True), Literal(Q.real, False))),
+        frozenset((Literal(Q.negative, True), Literal(Q.zero, True))),
+        frozenset((Literal(Q.nonnegative, False), Literal(Q.positive, True))),
+        frozenset((Literal(Q.nonnegative, False), Literal(Q.zero, True))),
+        frozenset((Literal(Q.nonnegative, True), Literal(Q.positive, False), Literal(Q.zero, False))),
+        frozenset((Literal(Q.nonpositive, False), Literal(Q.zero, True))),
+        frozenset((Literal(Q.nonzero, False), Literal(Q.positive, True))),
+        frozenset((Literal(Q.normal, False), Literal(Q.unitary, True))),
+        frozenset((Literal(Q.normal, True), Literal(Q.square, False))),
+        frozenset((Literal(Q.orthogonal, False), Literal(Q.real, True), Literal(Q.unitary, True))),
+        frozenset((Literal(Q.orthogonal, True), Literal(Q.positive_definite, False))),
+        frozenset((Literal(Q.orthogonal, True), Literal(Q.unitary, False))),
+        frozenset((Literal(Q.positive, False), Literal(Q.prime, True))),
+        frozenset((Literal(Q.positive, True), Literal(Q.real, False))),
+        frozenset((Literal(Q.positive, True), Literal(Q.zero, True))),
+        frozenset((Literal(Q.rational, True), Literal(Q.real, False))),
+        frozenset((Literal(Q.real, False), Literal(Q.zero, True))),
+        frozenset((Literal(Q.square, False), Literal(Q.symmetric, True))),
+        frozenset((Literal(Q.triangular, False), Literal(Q.unit_triangular, True))),
+        frozenset((Literal(Q.triangular, False), Literal(Q.upper_triangular, True)))
+    }
+
 # -{ Known facts in Conjunctive Normal Form }-
 @cacheit
 def get_known_facts_cnf():
@@ -85,7 +164,6 @@ def get_known_facts_cnf():
         Q.orthogonal | ~Q.real | ~Q.unitary,
         Q.negative | Q.positive | Q.zero | ~Q.real,
         Q.algebraic | Q.transcendental | ~Q.complex | ~Q.finite,
-        Q.composite | Q.prime | ~Q.integer | ~Q.positive,
         Q.irrational | Q.rational | ~Q.finite | ~Q.real
     )
 

sympy/assumptions/cnf.py

@@ -0,0 +1,378 @@
+"""
+The classes used here are for the internal use of assumptions system
+only and should not be used anywhere else as these don't possess the
+signatures common to SymPy objects. For general use of logic constructs
+please refer to sympy.logic classes And, Or, Not, etc.
+"""
+from itertools import combinations, product
+from sympy import S, Nor, Nand, Xor, Implies, Equivalent, ITE
+from sympy.logic.boolalg import Or, And, Not, Xnor
+from sympy.core.compatibility import zip_longest
+
+
+class Literal(object):
+    """
+    The smallest element of a CNF object
+    """
+
+    def __new__(cls, lit, is_Not=False):
+        if isinstance(lit, Not):
+            lit = lit.args[0]
+            is_Not = True
+        elif isinstance(lit, (AND, OR, Literal)):
+            return ~lit if is_Not else lit
+        obj = super(Literal, cls).__new__(cls)
+        obj.lit = lit
+        obj.is_Not = is_Not
+        return obj
+
+    @property
+    def arg(self):
+        return self.lit
+
+    def rcall(self, expr):
+        if callable(self.lit):
+            lit = self.lit(expr)
+        else:
+            try:
+                lit = self.lit.apply(expr)
+            except AttributeError:
+                lit = self.lit.rcall(expr)
+        return type(self)(lit, self.is_Not)
+
+    def __invert__(self):
+        is_Not = not self.is_Not
+        return Literal(self.lit, is_Not)
+
+    def __str__(self):
+        return '%s(%s, %s)' % (type(self).__name__, self.lit, self.is_Not)
+
+    __repr__ = __str__
+
+    def __eq__(self, other):
+        return self.arg == other.arg and self.is_Not == other.is_Not
+
+    def __hash__(self):
+        h = hash((type(self).__name__, self.arg, self.is_Not))
+        return h
+
+
+class OR(object):
+    """
+    A low-level implementation for Or
+    """
+    def __init__(self, *args):
+        self._args = args
+
+    @property
+    def args(self):
+        return sorted(self._args, key=str)
+
+    def rcall(self, expr):
+        return type(self)(*[arg.rcall(expr)
+                            for arg in self._args
+                            ])
+
+    def __invert__(self):
+        return AND(*[~arg for arg in self._args])
+
+    def __hash__(self):
+        return hash((type(self).__name__,) + tuple(self.args))
+
+    def __eq__(self, other):
+        return self.args == other.args
+
+    def __str__(self):
+        s = '(' + ' | '.join([str(arg) for arg in self.args]) + ')'
+        return s
+
+    __repr__ = __str__
+
+
+class AND(object):
+    """
+    A low-level implementation for And
+    """
+    def __init__(self, *args):
+        self._args = args
+
+    def __invert__(self):
+        return OR(*[~arg for arg in self._args])
+
+    @property
+    def args(self):
+        return sorted(self._args, key=str)
+
+    def rcall(self, expr):
+        return type(self)(*[arg.rcall(expr)
+                            for arg in self._args
+                            ])
+
+    def __hash__(self):
+        return hash((type(self).__name__,) + tuple(self.args))
+
+    def __eq__(self, other):
+        return self.args == other.args
+
+    def __str__(self):
+        s = '('+' & '.join([str(arg) for arg in self.args])+')'
+        return s
+
+    __repr__ = __str__
+
+
+def to_NNF(expr):
+    """
+    Generates the Negation Normal Form of any boolean expression in terms
+    of AND, OR, and Literal objects.
+    """
+
+    if isinstance(expr, Not):
+        arg = expr.args[0]
+        tmp = to_NNF(arg)  # Strategy: negate the NNF of expr
+        return ~tmp
+
+    if isinstance(expr, Or):
+        return OR(*[to_NNF(x) for x in Or.make_args(expr)])
+
+    if isinstance(expr, And):
+        return AND(*[to_NNF(x) for x in And.make_args(expr)])
+
+    if isinstance(expr, Nand):
+        tmp = AND(*[to_NNF(x) for x in expr.args])
+        return ~tmp
+
+    if isinstance(expr, Nor):
+        tmp = OR(*[to_NNF(x) for x in expr.args])
+        return ~tmp
+
+    if isinstance(expr, Xor):
+        cnfs = []
+        for i in range(0, len(expr.args) + 1, 2):
+            for neg in combinations(expr.args, i):
+                clause = [~to_NNF(s) if s in neg else to_NNF(s)
+                          for s in expr.args]
+                cnfs.append(OR(*clause))
+        return AND(*cnfs)
+
+    if isinstance(expr, Xnor):
+        cnfs = []
+        for i in range(0, len(expr.args) + 1, 2):
+            for neg in combinations(expr.args, i):
+                clause = [~to_NNF(s) if s in neg else to_NNF(s)
+                          for s in expr.args]
+                cnfs.append(OR(*clause))
+        return ~AND(*cnfs)
+
+    if isinstance(expr, Implies):
+        L, R = to_NNF(expr.args[0]), to_NNF(expr.args[1])
+        return OR(~L, R)
+
+    if isinstance(expr, Equivalent):
+        cnfs = []
+        for a, b in zip_longest(expr.args, expr.args[1:], fillvalue=expr.args[0]):
+            a = to_NNF(a)
+            b = to_NNF(b)
+            cnfs.append(OR(~a, b))
+        return AND(*cnfs)
+
+    if isinstance(expr, ITE):
+        L = to_NNF(expr.args[0])
+        M = to_NNF(expr.args[1])
+        R = to_NNF(expr.args[2])
+        return AND(OR(~L, M), OR(L, R))
+
+    else:
+        return Literal(expr)
+
+
+def distribute_AND_over_OR(expr):
+    """
+    Distributes AND over OR in the NNF expression.
+    Returns the result( Conjunctive Normal Form of expression)
+    as a CNF object.
+    """
+    if not isinstance(expr, (AND, OR)):
+        tmp = set()
+        tmp.add(frozenset((expr,)))
+        return CNF(tmp)
+
+    if isinstance(expr, OR):
+        return CNF.all_or(*[distribute_AND_over_OR(arg)
+                            for arg in expr._args])
+
+    if isinstance(expr, AND):
+        return CNF.all_and(*[distribute_AND_over_OR(arg)
+                             for arg in expr._args])
+
+
+class CNF(object):
+    """
+    Class to represent CNF of a Boolean expression.
+    Consists of set of clauses, which themselves are stored as
+    frozenset of Literal objects.
+    """
+    def __init__(self, clauses=None):
+        if not clauses:
+            clauses = set()
+        self.clauses = clauses
+
+    def add(self, prop):
+        clauses = CNF.to_CNF(prop).clauses
+        self.add_clauses(clauses)
+
+    def __str__(self):
+        s = ' & '.join(
+            ['(' + ' | '.join([str(lit) for lit in clause]) +')'
+            for clause in self.clauses]
+        )
+        return s
+
+    def extend(self, props):
+        for p in props:
+            self.add(p)
+        return self
+
+    def copy(self):
+        return CNF(set(self.clauses))
+
+    def add_clauses(self, clauses):
+        self.clauses |= clauses
+
+    @classmethod
+    def from_prop(cls, prop):
+        res = cls()
+        res.add(prop)
+        return res
+
+    def __iand__(self, other):
+        self.add_clauses(other.clauses)
+        return self
+
+    def all_predicates(self):
+        predicates = set()
+        for c in self.clauses:
+            predicates |= {arg.lit for arg in c}
+        return predicates
+
+    def _or(self, cnf):
+        clauses = set()
+        for a, b in product(self.clauses, cnf.clauses):
+            tmp = set(a)
+            for t in b:
+                tmp.add(t)
+            clauses.add(frozenset(tmp))
+        return CNF(clauses)
+
+    def _and(self, cnf):
+        clauses = self.clauses.union(cnf.clauses)
+        return CNF(clauses)
+
+    def _not(self):
+        clss = list(self.clauses)
+        ll = set()
+        for x in clss[-1]:
+            ll.add(frozenset((~x,)))
+        ll = CNF(ll)
+
+        for rest in clss[:-1]:
+            p = set()
+            for x in rest:
+                p.add(frozenset((~x,)))
+            ll = ll._or(CNF(p))
+        return ll
+
+    def rcall(self, expr):
+        clause_list = list()
+        for clause in self.clauses:
+            lits = [arg.rcall(expr) for arg in clause]
+            clause_list.append(OR(*lits))
+        expr = AND(*clause_list)
+        return distribute_AND_over_OR(expr)
+
+
+
+    @classmethod
+    def all_or(cls, *cnfs):
+        b = cnfs[0].copy()
+        for rest in cnfs[1:]:
+            b = b._or(rest)
+        return b
+
+    @classmethod
+    def all_and(cls, *cnfs):
+        b = cnfs[0].copy()
+        for rest in cnfs[1:]:
+            b = b._and(rest)
+        return b
+
+    @classmethod
+    def to_CNF(cls, expr):
+        expr = to_NNF(expr)
+        expr = distribute_AND_over_OR(expr)
+        return expr
+
+    @classmethod
+    def CNF_to_cnf(cls, cnf):
+        """
+        Converts CNF object to SymPy's boolean expression
+        retaining the form of expression.
+        """
+        def remove_literal(arg):
+            return Not(arg.lit) if arg.is_Not else arg.lit
+
+        return And(*(Or(*(remove_literal(arg) for arg in clause)) for clause in cnf.clauses))
+
+
+class EncodedCNF(object):
+    """
+    Class for encoding the CNF expression.
+    """
+    def __init__(self, data=None, encoding=None):
+        if not data and not encoding:
+            data = list()
+            encoding = dict()
+        self.data = data
+        self.encoding = encoding
+        self._symbols = list(encoding.keys())
+
+    def from_cnf(self, cnf):
+        self._symbols = list(cnf.all_predicates())
+        n = len(self._symbols)
+        self.encoding = dict(list(zip(self._symbols, list(range(1, n + 1)))))
+        self.data = [self.encode(clause) for clause in cnf.clauses]
+
+    @property
+    def symbols(self):
+        return self._symbols
+
+    @property
+    def variables(self):
+        return range(1, len(self._symbols) + 1)
+
+    def copy(self):
+        new_data = [set(clause) for clause in self.data]
+        return EncodedCNF(new_data, dict(self.encoding))
+
+    def add_prop(self, prop):
+        cnf = CNF.from_prop(prop)
+        self.add_from_cnf(cnf)
+
+    def add_from_cnf(self, cnf):
+        clauses = [self.encode(clause) for clause in cnf.clauses]
+        self.data += clauses
+
+    def encode_arg(self, arg):
+        literal = arg.lit
+        value = self.encoding.get(literal, None)
+        if value is None:
+            n = len(self._symbols)
+            self._symbols.append(literal)
+            value = self.encoding[literal] = n + 1
+        if arg.is_Not:
+            return -value
+        else:
+            return value
+
+    def encode(self, clause):
+        return {self.encode_arg(arg) if not arg.lit == S.false else 0 for arg in clause}

sympy/assumptions/handlers/sets.py

@@ -6,8 +6,9 @@ from __future__ import print_function, division
 from sympy.assumptions import Q, ask
 from sympy.assumptions.handlers import CommonHandler, test_closed_group
 from sympy.core.numbers import pi
+from sympy.core.logic import fuzzy_bool
 from sympy.functions.elementary.exponential import exp, log
-from sympy import I
+from sympy import I, Eq, conjugate, MatrixBase
 
 
 class AskIntegerHandler(CommonHandler):
@@ -336,6 +337,12 @@ class AskHermitianHandler(AskRealHandler):
     Test that an expression belongs to the field of Hermitian operators
     """
 
+    @staticmethod
+    def Expr(expr, assumptions):
+        if isinstance(expr, MatrixBase):
+            return None
+        return AskRealHandler.Expr(expr, assumptions)
+
     @staticmethod
     def Add(expr, assumptions):
         """
@@ -388,6 +395,18 @@ class AskHermitianHandler(AskRealHandler):
 
     cos, exp = [sin]*2
 
+    @staticmethod
+    def MatrixBase(mat, assumptions):
+        rows, cols = mat.shape
+        ret_val = True
+        for i in range(rows):
+            for j in range(i, cols):
+                cond = fuzzy_bool(Eq(mat[i, j], conjugate(mat[j, i])))
+                if cond == None:
+                    ret_val = None
+                if cond == False:
+                    return False
+        return ret_val
 
 class AskComplexHandler(CommonHandler):
     """

sympy/assumptions/refine.py

@@ -240,6 +240,57 @@ def refine_Relational(expr, assumptions):
     return ask(Q.is_true(expr), assumptions)
 
 
+def refine_re(expr, assumptions):
+    """
+    Handler for real part.
+
+    >>> from sympy.assumptions.refine import refine_re
+    >>> from sympy import Q, re
+    >>> from sympy.abc import x
+    >>> refine_re(re(x), Q.real(x))
+    x
+    >>> refine_re(re(x), Q.imaginary(x))
+    0
+    """
+    arg = expr.args[0]
+    if ask(Q.real(arg), assumptions):
+        return arg
+    if ask(Q.imaginary(arg), assumptions):
+        return 0
+    return _refine_reim(expr, assumptions)
+
+
+def refine_im(expr, assumptions):
+    """
+    Handler for imaginary part.
+
+    >>> from sympy.assumptions.refine import refine_im
+    >>> from sympy import Q, im
+    >>> from sympy.abc import x
+    >>> refine_im(im(x), Q.real(x))
+    0
+    >>> refine_im(im(x), Q.imaginary(x))
+    -I*x
+    """
+    arg = expr.args[0]
+    if ask(Q.real(arg), assumptions):
+        return 0
+    if ask(Q.imaginary(arg), assumptions):
+        return - S.ImaginaryUnit * arg
+    return _refine_reim(expr, assumptions)
+
+
+def _refine_reim(expr, assumptions):
+    # Helper function for refine_re & refine_im
+    expanded = expr.expand(complex = True)
+    if expanded != expr:
+        refined = refine(expanded, assumptions)
+        if refined != expanded:
+            return refined
+    # Best to leave the expression as is
+    return None
+
+
 handlers_dict = {
     'Abs': refine_abs,
     'Pow': refine_Pow,
@@ -249,5 +300,7 @@ handlers_dict = {
     'GreaterThan': refine_Relational,
     'LessThan': refine_Relational,
     'StrictGreaterThan': refine_Relational,
-    'StrictLessThan': refine_Relational
+    'StrictLessThan': refine_Relational,
+    're': refine_re,
+    'im': refine_im
 }

sympy/assumptions/satask.py

@@ -1,22 +1,39 @@
 from __future__ import print_function, division
 
-from sympy.assumptions.ask_generated import get_known_facts_cnf
+from sympy import Symbol, S
+from sympy.assumptions.ask_generated import get_all_known_facts
 from sympy.assumptions.assume import global_assumptions, AppliedPredicate
 from sympy.assumptions.sathandlers import fact_registry
-from sympy.core import oo, Tuple
+from sympy.core import oo
 from sympy.logic.inference import satisfiable
-from sympy.logic.boolalg import And
+from sympy.assumptions.cnf import CNF, EncodedCNF
 
 
 def satask(proposition, assumptions=True, context=global_assumptions,
-    use_known_facts=True, iterations=oo):
-    relevant_facts = get_all_relevant_facts(proposition, assumptions, context,
-        use_known_facts=use_known_facts, iterations=iterations)
+        use_known_facts=True, iterations=oo):
+    props = CNF.from_prop(proposition)
+    _props = CNF.from_prop(~proposition)
+    if context:
+        tmp = CNF()
+        context = tmp.extend(context)
+    assumptions = CNF.from_prop(assumptions)
 
-    can_be_true = satisfiable(And(proposition, assumptions,
-        relevant_facts, *context))
-    can_be_false = satisfiable(And(~proposition, assumptions,
-        relevant_facts, *context))
+    sat = get_all_relevant_facts(props, assumptions, context,
+        use_known_facts=use_known_facts, iterations=iterations)
+    if context:
+        sat.add_from_cnf(context)
+    sat.add_from_cnf(assumptions)
+
+    return check_satisfiability(props, _props, sat)
+
+
+def check_satisfiability(prop, _prop, factbase):
+    sat_true = factbase.copy()
+    sat_false = factbase.copy()
+    sat_true.add_from_cnf(prop)
+    sat_false.add_from_cnf(_prop)
+    can_be_true = satisfiable(sat_true)
+    can_be_false = satisfiable(sat_false)
 
     if can_be_true and can_be_false:
         return None
@@ -34,35 +51,61 @@ def satask(proposition, assumptions=True, context=global_assumptions,
         raise ValueError("Inconsistent assumptions")
 
 
-def get_relevant_facts(proposition, assumptions=(True,),
-    context=global_assumptions, use_known_facts=True, exprs=None,
+def get_relevant_facts(proposition, assumptions=None,
+    context=None, exprs=None,
     relevant_facts=None):
 
     newexprs = set()
-    if not exprs:
-        keys = proposition.atoms(AppliedPredicate)
-        # XXX: We need this since True/False are not Basic
-        keys |= Tuple(*assumptions).atoms(AppliedPredicate)
-        if context:
-            keys |= And(*context).atoms(AppliedPredicate)
 
-        exprs = {key.args[0] for key in keys}
+    if not assumptions:
+        assumptions = CNF({S.true})
 
     if not relevant_facts:
-        relevant_facts = set([])
+        relevant_facts = set()
 
-    if use_known_facts:
-        for expr in exprs:
-            relevant_facts.add(get_known_facts_cnf().rcall(expr))
+    def find_symbols(pred):
+        if isinstance(pred, CNF):
+            symbols = set()
+            for a in pred.all_predicates():
+                symbols |= find_symbols(a)
+            return symbols
+        if isinstance(pred.args, AppliedPredicate):
+            return {pred.args[0]}
+        return pred.atoms(Symbol)
+
+    if not exprs:
+        req_keys = find_symbols(proposition)
+        keys = proposition.all_predicates()
+        # XXX: We need this since True/False are not Basic
+        lkeys = set()
+        lkeys |= assumptions.all_predicates()
+        if context:
+            lkeys |= context.all_predicates()
+
+        lkeys = lkeys - {S.true, S.false}
+        tmp_keys = None
+        while tmp_keys != set():
+            tmp = set()
+            for l in lkeys:
+                syms = find_symbols(l)
+                if (syms & req_keys) != set():
+                    tmp |= syms
+            tmp_keys = tmp - req_keys
+            req_keys |= tmp_keys
+        keys |= {l for l in lkeys if find_symbols(l) & req_keys != set()}
+
+        exprs = {key.args[0] if isinstance(key, AppliedPredicate) else key for key in keys}
+        return exprs, relevant_facts
 
     for expr in exprs:
         for fact in fact_registry[expr.func]:
-            newfact = fact.rcall(expr)
-            relevant_facts.add(newfact)
-            newexprs |= set([key.args[0] for key in
-                newfact.atoms(AppliedPredicate)])
+            cnf_fact = CNF.to_CNF(fact)
+            newfact = cnf_fact.rcall(expr)
+            relevant_facts = relevant_facts._and(newfact)
+            newexprs |= set([key.args[0] for key in newfact.all_predicates()
+                             if isinstance(key, AppliedPredicate)])
 
-    return relevant_facts, newexprs - exprs
+    return newexprs - exprs, relevant_facts
 
 
 def get_all_relevant_facts(proposition, assumptions=True,
@@ -72,15 +115,44 @@ def get_all_relevant_facts(proposition, assumptions=True,
     # we stop getting new things. Hopefully this strategy won't lead to an
     # infinite loop in the future.
     i = 0
-    relevant_facts = set()
+    relevant_facts = CNF()
     exprs = None
+    all_exprs = set()
     while exprs != set():
-        (relevant_facts, exprs) = get_relevant_facts(proposition,
-                And.make_args(assumptions), context,
-                use_known_facts=use_known_facts, exprs=exprs,
+        exprs, relevant_facts = get_relevant_facts(proposition,
+                assumptions, context, exprs=exprs,
                 relevant_facts=relevant_facts)
+        all_exprs |= exprs
         i += 1
         if i >= iterations:
-            return And(*relevant_facts)
+            break
 
-    return And(*relevant_facts)
+    if use_known_facts:
+        known_facts_CNF = CNF()
+        known_facts_CNF.add_clauses(get_all_known_facts())
+        kf_encoded = EncodedCNF()
+        kf_encoded.from_cnf(known_facts_CNF)
+
+        def translate_literal(lit, delta):
+            if lit > 0:
+                return lit + delta
+            else:
+                return lit - delta
+
+        def translate_data(data, delta):
+            return [{translate_literal(i, delta) for i in clause} for clause in data]
+        data = []
+        symbols = []
+        n_lit = len(kf_encoded.symbols)
+        for i, expr in enumerate(all_exprs):
+            symbols += [pred(expr) for pred in kf_encoded.symbols]
+            data += translate_data(kf_encoded.data, i * n_lit)
+
+        encoding = dict(list(zip(symbols, range(1, len(symbols)+1))))
+        ctx = EncodedCNF(data, encoding)
+    else:
+        ctx = EncodedCNF()
+
+    ctx.add_from_cnf(relevant_facts)
+
+    return ctx

sympy/assumptions/sathandlers.py

@@ -4,20 +4,19 @@ from collections import defaultdict
 
 from sympy.assumptions.ask import Q
 from sympy.assumptions.assume import Predicate, AppliedPredicate
+from sympy.assumptions.cnf import AND, OR, to_NNF
 from sympy.core import (Add, Mul, Pow, Integer, Number, NumberSymbol,)
 from sympy.core.compatibility import MutableMapping
 from sympy.core.numbers import ImaginaryUnit
-from sympy.core.logic import fuzzy_or, fuzzy_and
 from sympy.core.rules import Transform
 from sympy.core.sympify import _sympify
 from sympy.functions.elementary.complexes import Abs
-from sympy.logic.boolalg import (Equivalent, Implies, And, Or,
-    BooleanFunction, Not)
+from sympy.logic.boolalg import (Equivalent, Implies, BooleanFunction)
 from sympy.matrices.expressions import MatMul
 
 # APIs here may be subject to change
 
-# XXX: Better name?
+
 class UnevaluatedOnFree(BooleanFunction):
     """
     Represents a Boolean function that remains unevaluated on free predicates
@@ -68,13 +67,19 @@ class UnevaluatedOnFree(BooleanFunction):
         obj.expr = predicate_args.pop()
         obj.pred = arg.xreplace(Transform(lambda e: e.func, lambda e:
             isinstance(e, AppliedPredicate)))
-        applied = obj.apply()
+        applied = obj.apply(obj.expr)
         if applied is None:
             return obj
         return applied
 
-    def apply(self):
-        return
+    def apply(self, expr=None):
+        if expr is None:
+            return
+        pred = to_NNF(self.pred)
+        return self._eval_apply(expr, pred)
+
+    def _eval_apply(self, expr, pred):
+        return None
 
 
 class AllArgs(UnevaluatedOnFree):
@@ -89,7 +94,6 @@ class AllArgs(UnevaluatedOnFree):
 
     Example
     =======
-
     >>> from sympy.assumptions.sathandlers import AllArgs
     >>> from sympy import symbols, Q
     >>> x, y = symbols('x y')
@@ -97,11 +101,13 @@ class AllArgs(UnevaluatedOnFree):
     >>> a
     AllArgs(Q.negative | Q.positive)
     >>> a.rcall(x*y)
-    (Q.negative(x) | Q.positive(x)) & (Q.negative(y) | Q.positive(y))
+    ((Literal(Q.negative(x), False) | Literal(Q.positive(x), False)) & (Literal(Q.negative(y), False) | \
+    Literal(Q.positive(y), False)))
+
     """
 
-    def apply(self):
-        return And(*[self.pred.rcall(arg) for arg in self.expr.args])
+    def _eval_apply(self, expr, pred):
+        return AND(*[pred.rcall(arg) for arg in expr.args])
 
 
 class AnyArgs(UnevaluatedOnFree):
@@ -116,7 +122,6 @@ class AnyArgs(UnevaluatedOnFree):
 
     Example
     =======
-
     >>> from sympy.assumptions.sathandlers import AnyArgs
     >>> from sympy import symbols, Q
     >>> x, y = symbols('x y')
@@ -124,11 +129,13 @@ class AnyArgs(UnevaluatedOnFree):
     >>> a
     AnyArgs(Q.negative & Q.positive)
     >>> a.rcall(x*y)
-    (Q.negative(x) & Q.positive(x)) | (Q.negative(y) & Q.positive(y))
+    ((Literal(Q.negative(x), False) & Literal(Q.positive(x), False)) | (Literal(Q.negative(y), False) & \
+    Literal(Q.positive(y), False)))
+
     """
 
-    def apply(self):
-        return Or(*[self.pred.rcall(arg) for arg in self.expr.args])
+    def _eval_apply(self, expr, pred):
+        return OR(*[pred.rcall(arg) for arg in expr.args])
 
 
 class ExactlyOneArg(UnevaluatedOnFree):
@@ -144,7 +151,6 @@ class ExactlyOneArg(UnevaluatedOnFree):
 
     Example
     =======
-
     >>> from sympy.assumptions.sathandlers import ExactlyOneArg
     >>> from sympy import symbols, Q
     >>> x, y = symbols('x y')
@@ -152,17 +158,19 @@ class ExactlyOneArg(UnevaluatedOnFree):
     >>> a
     ExactlyOneArg(Q.positive)
     >>> a.rcall(x*y)
-    (Q.positive(x) & ~Q.positive(y)) | (Q.positive(y) & ~Q.positive(x))
+    ((Literal(Q.positive(x), False) & Literal(Q.positive(y), True)) | (Literal(Q.positive(x), True) & \
+    Literal(Q.positive(y), False)))
+
     """
-    def apply(self):
-        expr = self.expr
-        pred = self.pred
+
+    def _eval_apply(self, expr, pred):
         pred_args = [pred.rcall(arg) for arg in expr.args]
         # Technically this is xor, but if one term in the disjunction is true,
         # it is not possible for the remainder to be true, so regular or is
         # fine in this case.
-        return Or(*[And(pred_args[i], *map(Not, pred_args[:i] +
-            pred_args[i+1:])) for i in range(len(pred_args))])
+        res = OR(*[AND(pred_args[i], *[~lit for lit in pred_args[:i] +
+            pred_args[i+1:]]) for i in range(len(pred_args))])
+        return res
         # Note: this is the equivalent cnf form. The above is more efficient
         # as the first argument of an implication, since p >> q is the same as
         # q | ~p, so the the ~ will convert the Or to and, and one just needs
@@ -172,10 +180,6 @@ class ExactlyOneArg(UnevaluatedOnFree):
 
 
 def _old_assump_replacer(obj):
-    # Things to be careful of:
-    # - real means real or infinite in the old assumptions.
-    # - nonzero does not imply real in the old assumptions.
-    # - finite means finite and not zero in the old assumptions.
     if not isinstance(obj, AppliedPredicate):
         return obj
 
@@ -183,33 +187,34 @@ def _old_assump_replacer(obj):
     ret = None
 
     if obj.func == Q.positive:
-        ret = fuzzy_and([e.is_finite, e.is_positive])
-    if obj.func == Q.zero:
+        ret = e.is_positive
+    elif obj.func == Q.zero:
         ret = e.is_zero
-    if obj.func == Q.negative:
-        ret = fuzzy_and([e.is_finite, e.is_negative])
-    if obj.func == Q.nonpositive:
-        ret = fuzzy_and([e.is_finite, e.is_nonpositive])
-    if obj.func == Q.nonzero:
-        ret = fuzzy_and([e.is_nonzero, e.is_finite])
-    if obj.func == Q.nonnegative:
-        ret = fuzzy_and([fuzzy_or([e.is_zero, e.is_finite]),
-        e.is_nonnegative])
+    elif obj.func == Q.negative:
+        ret = e.is_negative
+    elif obj.func == Q.nonpositive:
+        ret = e.is_nonpositive
+    elif obj.func == Q.nonzero:
+        ret = e.is_nonzero
+    elif obj.func == Q.nonnegative:
+        ret = e.is_nonnegative
 
-    if obj.func == Q.rational:
+    elif obj.func == Q.rational:
         ret = e.is_rational
-    if obj.func == Q.irrational:
+    elif obj.func == Q.irrational:
         ret = e.is_irrational
 
-    if obj.func == Q.even:
+    elif obj.func == Q.even:
         ret = e.is_even
-    if obj.func == Q.odd:
+    elif obj.func == Q.odd:
         ret = e.is_odd
-    if obj.func == Q.integer:
+    elif obj.func == Q.integer:
         ret = e.is_integer
-    if obj.func == Q.imaginary:
+    elif obj.func == Q.composite:
+        ret = e.is_composite
+    elif obj.func == Q.imaginary:
         ret = e.is_imaginary
-    if obj.func == Q.commutative:
+    elif obj.func == Q.commutative:
         ret = e.is_commutative
 
     if ret is None:
@@ -229,15 +234,18 @@ def evaluate_old_assump(pred):
 
 
 class CheckOldAssump(UnevaluatedOnFree):
-    def apply(self):
-        return Equivalent(self.args[0], evaluate_old_assump(self.args[0]))
+    def apply(self, expr=None, is_Not=False):
+        arg = self.args[0](expr) if callable(self.args[0]) else self.args[0]
+        res = Equivalent(arg, evaluate_old_assump(arg))
+        return to_NNF(res)
 
 
 class CheckIsPrime(UnevaluatedOnFree):
-    def apply(self):
+    def apply(self, expr=None, is_Not=False):
         from sympy import isprime
-        return Equivalent(self.args[0], isprime(self.expr))
-
+        arg = self.args[0](expr) if callable(self.args[0]) else self.args[0]
+        res = Equivalent(arg, isprime(expr))
+        return to_NNF(res)
 
 class CustomLambda(object):
     """
@@ -248,8 +256,8 @@ class CustomLambda(object):
     def __init__(self, lamda):
         self.lamda = lamda
 
-    def rcall(self, *args):
-        return self.lamda(*args)
+    def apply(self, *args):
+        return to_NNF(self.lamda(*args))
 
 
 class ClassFactRegistry(MutableMapping):
@@ -313,6 +321,7 @@ for klass, fact in [
     # matching, so that we can just write Equivalent(Q.zero(x**y), Q.zero(x) & Q.positive(y))
     (Pow, CustomLambda(lambda power: Equivalent(Q.zero(power), Q.zero(power.base) & Q.positive(power.exp)))),
     (Integer, CheckIsPrime(Q.prime)),
+    (Integer, CheckOldAssump(Q.composite)),
     # Implicitly assumes Mul has more than one arg
     # Would be AllArgs(Q.prime | Q.composite) except 1 is composite
     (Mul, Implies(AllArgs(Q.prime), ~Q.prime)),

sympy/assumptions/tests/test_assumptions_2.py

@@ -5,9 +5,7 @@ from sympy.abc import x, y
 from sympy.assumptions.assume import global_assumptions, Predicate
 from sympy.assumptions.ask import _extract_facts, Q
 from sympy.core import symbols
-from sympy.logic.boolalg import Or
 from sympy.printing import pretty
-from sympy.utilities.pytest import XFAIL
 
 
 def test_equal():

sympy/assumptions/tests/test_query.py

@@ -16,10 +16,12 @@ from sympy.functions.elementary.miscellaneous import sqrt
 from sympy.functions.elementary.trigonometric import (
     acos, acot, asin, atan, cos, cot, sin, tan)
 from sympy.logic.boolalg import Equivalent, Implies, Xor, And, to_cnf
+from sympy.matrices import Matrix, SparseMatrix
 from sympy.utilities.pytest import XFAIL, slow, raises, warns_deprecated_sympy
 from sympy.assumptions.assume import assuming
 import math
 
+
 def test_int_1():
     z = 1
     assert ask(Q.commutative(z)) is True
@@ -174,7 +176,8 @@ def test_infinity():
     assert ask(Q.complex(oo)) is False
     assert ask(Q.irrational(oo)) is False
     assert ask(Q.imaginary(oo)) is False
-    assert ask(Q.positive(oo)) is True
+    assert ask(Q.positive(oo)) is False
+    #assert ask(Q.extended_positive(oo)) is True
     assert ask(Q.negative(oo)) is False
     assert ask(Q.even(oo)) is False
     assert ask(Q.odd(oo)) is False
@@ -197,7 +200,8 @@ def test_neg_infinity():
     assert ask(Q.irrational(mm)) is False
     assert ask(Q.imaginary(mm)) is False
     assert ask(Q.positive(mm)) is False
-    assert ask(Q.negative(mm)) is True
+    assert ask(Q.negative(mm)) is False
+    #assert ask(Q.extended_negative(mm)) is True
     assert ask(Q.even(mm)) is False
     assert ask(Q.odd(mm)) is False
     assert ask(Q.finite(mm)) is False
@@ -565,8 +569,8 @@ def test_I():
     assert ask(Q.real(z)) is True
 
 
-@slow
-def test_bounded1():
+
+def test_bounded():
     x, y, z = symbols('x,y,z')
     assert ask(Q.finite(x)) is None
     assert ask(Q.finite(x), Q.finite(x)) is True
@@ -652,9 +656,6 @@ def test_bounded1():
     assert ask(Q.finite(a), ~Q.positive(x) & Q.positive(y)) is None
     assert ask(Q.finite(a), ~Q.positive(x) & ~Q.positive(y)) is None
 
-
-@slow
-def test_bounded2a():
     x, y, z = symbols('x,y,z')
     a = x + y + z
     x, y, z = a.args
@@ -799,13 +800,6 @@ def test_bounded2a():
         Q.finite(a), Q.finite(x) & ~Q.finite(y) & Q.positive(z)) is None
     assert ask(Q.finite(a), Q.finite(x) & Q.positive(y) &
         ~Q.finite(y) & Q.positive(z) & ~Q.finite(z)) is False
-
-
-@slow
-def test_bounded2b():
-    x, y, z = symbols('x,y,z')
-    a = x + y + z
-    x, y, z = a.args
     assert ask(Q.finite(a), Q.finite(x) &
         Q.positive(y) & ~Q.finite(y) & Q.negative(z)) is None
     assert ask(
@@ -980,9 +974,6 @@ def test_bounded2b():
     assert ask(Q.finite(2*x)) is None
     assert ask(Q.finite(2*x), Q.finite(x)) is True
 
-
-@slow
-def test_bounded3():
     x, y, z = symbols('x,y,z')
     a = x*y
     x, y = a.args
@@ -1041,9 +1032,9 @@ def test_bounded3():
     assert ask(Q.finite(2**x)) is None
     assert ask(Q.finite(2**x), Q.finite(x)) is True
     assert ask(Q.finite(x**x)) is None
-    assert ask(Q.finite(Rational(1, 2) ** x)) is None
-    assert ask(Q.finite(Rational(1, 2) ** x), Q.positive(x)) is True
-    assert ask(Q.finite(Rational(1, 2) ** x), Q.negative(x)) is None
+    assert ask(Q.finite(S.Half ** x)) is None
+    assert ask(Q.finite(S.Half ** x), Q.positive(x)) is True
+    assert ask(Q.finite(S.Half ** x), Q.negative(x)) is None
     assert ask(Q.finite(2**x), Q.negative(x)) is True
     assert ask(Q.finite(sqrt(x))) is None
     assert ask(Q.finite(2**x), ~Q.finite(x)) is False
@@ -1176,7 +1167,6 @@ def test_complex():
     assert ask(Q.complex(im(x))) is True
 
 
-@slow
 def test_even_query():
     assert ask(Q.even(x)) is None
     assert ask(Q.even(x), Q.integer(x)) is None
@@ -1331,7 +1321,6 @@ def test_rational():
         assert ask(Q.rational(h(x)), Q.rational(x)) is False
 
 
-@slow
 def test_hermitian():
     assert ask(Q.hermitian(x)) is None
     assert ask(Q.hermitian(x), Q.antihermitian(x)) is False
@@ -1435,7 +1424,6 @@ def test_hermitian():
         Q.imaginary(x) & Q.imaginary(y) & Q.imaginary(z)) is True
 
 
-@slow
 def test_imaginary():
     assert ask(Q.imaginary(x)) is None
     assert ask(Q.imaginary(x), Q.real(x)) is False
@@ -1512,7 +1500,7 @@ def test_imaginary():
     assert ask(Q.imaginary(exp(pi*I/2, evaluate=False))) is True
 
     # issue 7886
-    assert ask(Q.imaginary(Pow(x, S.One/4)), Q.real(x) & Q.negative(x)) is False
+    assert ask(Q.imaginary(Pow(x, Rational(1, 4))), Q.real(x) & Q.negative(x)) is False
 
 
 def test_integer():
@@ -1537,7 +1525,6 @@ def test_integer():
     assert ask(Q.integer(x/3), Q.even(x)) is None
 
 
-@slow
 def test_negative():
     assert ask(Q.negative(x), Q.negative(x)) is True
     assert ask(Q.negative(x), Q.positive(x)) is False
@@ -1609,7 +1596,6 @@ def test_nonzero():
     assert ask(Q.nonzero(cos(1)**2 + sin(1)**2 - 1)) is None
 
 
-@slow
 def test_zero():
     assert ask(Q.zero(x)) is None
     assert ask(Q.zero(x), Q.real(x)) is None
@@ -1640,7 +1626,6 @@ def test_zero():
     assert ask(Q.zero(x) | Q.zero(y), Q.zero(x*y)) is True
 
 
-@slow
 def test_odd_query():
     assert ask(Q.odd(x)) is None
     assert ask(Q.odd(x), Q.odd(x)) is True
@@ -1746,7 +1731,6 @@ def test_prime():
     assert ask(Q.prime(x**y), Q.integer(x) & Q.integer(y)) is False
 
 
-@slow
 def test_positive():
     assert ask(Q.positive(x), Q.positive(x)) is True
     assert ask(Q.positive(x), Q.negative(x)) is False
@@ -1846,7 +1830,7 @@ def test_real_basic():
     assert ask(Q.real(I*x), Q.imaginary(x)) is True
     assert ask(Q.real(I*x), Q.complex(x)) is None
 
-@slow
+
 def test_real_pow():
     assert ask(Q.real(x**2), Q.real(x)) is True
     assert ask(Q.real(sqrt(x)), Q.negative(x)) is False
@@ -1874,6 +1858,7 @@ def test_real_pow():
     assert ask(Q.real(x**i), Q.imaginary(i)) is None  # x could be 0
     assert ask(Q.real(x**(I*pi/log(x))), Q.real(x)) is True
 
+
 def test_real_functions():
     # trigonometric functions
     assert ask(Q.real(sin(x))) is None
@@ -1907,6 +1892,18 @@ def test_real_functions():
     assert ask(Q.real(im(x))) is True
 
 
+def test_matrix():
+
+    # hermitian
+    assert ask(Q.hermitian(Matrix([[2, 2 + I, 4], [2 - I, 3, I], [4, -I, 1]]))) == True
+    assert ask(Q.hermitian(Matrix([[2, 2 + I, 4], [2 + I, 3, I], [4, -I, 1]]))) == False
+    z = symbols('z', complex=True)
+    assert ask(Q.hermitian(Matrix([[2, 2 + I, z], [2 - I, 3, I], [4, -I, 1]]))) == None
+    assert ask(Q.hermitian(SparseMatrix(((25, 15, -5), (15, 18, 0), (-5, 0, 11))))) == True
+    assert ask(Q.hermitian(SparseMatrix(((25, 15, -5), (15, I, 0), (-5, 0, 11))))) == False
+    assert ask(Q.hermitian(SparseMatrix(((25, 15, -5), (15, z, 0), (-5, 0, 11))))) == None
+
+
 def test_algebraic():
     assert ask(Q.algebraic(x)) is None
 
@@ -1921,8 +1918,8 @@ def test_algebraic():
     assert ask(Q.algebraic(I*sqrt(3))) is True
     assert ask(Q.algebraic(sqrt(1 + I*sqrt(3)))) is True
 
-    assert ask(Q.algebraic((1 + I*sqrt(3)**(S(17)/31)))) is True
-    assert ask(Q.algebraic((1 + I*sqrt(3)**(S(17)/pi)))) is False
+    assert ask(Q.algebraic((1 + I*sqrt(3)**Rational(17, 31)))) is True
+    assert ask(Q.algebraic((1 + I*sqrt(3)**(17/pi)))) is False
 
     for f in [exp, sin, tan, asin, atan, cos]:
         assert ask(Q.algebraic(f(7))) is False
@@ -2009,7 +2006,6 @@ def test_composite_assumptions():
     assert ask(Q.positive(x), Q.real(x) >> Q.positive(y)) is None
     assert ask(Q.real(x), ~(Q.real(x) >> Q.real(y))) is True
 
-
 def test_incompatible_resolutors():
     class Prime2AskHandler(AskHandler):
         @staticmethod
@@ -2182,7 +2178,7 @@ def test_check_old_assumption():
     assert ask(Q.real(x)) is None
     assert ask(Q.complex(x)) is True
 
-    x = symbols('x', positive=True)
+    x = symbols('x', positive=True, finite=True)
     assert ask(Q.positive(x)) is True
     assert ask(Q.negative(x)) is False
     assert ask(Q.real(x)) is True

sympy/assumptions/tests/test_refine.py

@@ -1,6 +1,6 @@
 from sympy import (Abs, exp, Expr, I, pi, Q, Rational, refine, S, sqrt,
-                   atan, atan2, nan, Symbol)
-from sympy.abc import x, y, z
+                   atan, atan2, nan, Symbol, re, im)
+from sympy.abc import w, x, y, z
 from sympy.core.relational import Eq, Ne
 from sympy.functions.elementary.piecewise import Piecewise
 from sympy.utilities.pytest import slow
@@ -26,17 +26,14 @@ def test_pow1():
     assert refine(sqrt(x**2), Q.complex(x)) != Abs(x)
     assert refine(sqrt(x**2), Q.real(x)) == Abs(x)
     assert refine(sqrt(x**2), Q.positive(x)) == x
-    assert refine((x**3)**(S(1)/3)) != x
+    assert refine((x**3)**Rational(1, 3)) != x
 
-    assert refine((x**3)**(S(1)/3), Q.real(x)) != x
-    assert refine((x**3)**(S(1)/3), Q.positive(x)) == x
+    assert refine((x**3)**Rational(1, 3), Q.real(x)) != x
+    assert refine((x**3)**Rational(1, 3), Q.positive(x)) == x
 
     assert refine(sqrt(1/x), Q.real(x)) != 1/sqrt(x)
     assert refine(sqrt(1/x), Q.positive(x)) == 1/sqrt(x)
 
-
-@slow
-def test_pow2():
     # powers of (-1)
     assert refine((-1)**(x + y), Q.even(x)) == (-1)**y
     assert refine((-1)**(x + y + z), Q.odd(x) & Q.odd(z)) == (-1)**y
@@ -44,9 +41,6 @@ def test_pow2():
     assert refine((-1)**(x + y + 2), Q.odd(x)) == (-1)**(y + 1)
     assert refine((-1)**(x + 3)) == (-1)**(x + 1)
 
-
-@slow
-def test_pow3():
     # continuation
     assert refine((-1)**((-1)**x/2 - S.Half), Q.integer(x)) == (-1)**x
     assert refine((-1)**((-1)**x/2 + S.Half), Q.integer(x)) == (-1)**(x + 1)
@@ -54,7 +48,7 @@ def test_pow3():
 
 
 @slow
-def test_pow4():
+def test_pow2():
     assert refine((-1)**((-1)**x/2 - 7*S.Half), Q.integer(x)) == (-1)**(x + 1)
     assert refine((-1)**((-1)**x/2 - 9*S.Half), Q.integer(x)) == (-1)**x
 
@@ -67,7 +61,7 @@ def test_pow4():
 def test_exp():
     x = Symbol('x', integer=True)
     assert refine(exp(pi*I*2*x)) == 1
-    assert refine(exp(pi*I*2*(x + Rational(1, 2)))) == -1
+    assert refine(exp(pi*I*2*(x + S.Half))) == -1
     assert refine(exp(pi*I*2*(x + Rational(1, 4)))) == I
     assert refine(exp(pi*I*2*(x + Rational(3, 4)))) == -I
 
@@ -142,7 +136,40 @@ def test_atan2():
     assert refine(atan2(y, x), Q.zero(y) & Q.negative(x)) == pi
     assert refine(atan2(y, x), Q.positive(y) & Q.zero(x)) == pi/2
     assert refine(atan2(y, x), Q.negative(y) & Q.zero(x)) == -pi/2
-    assert refine(atan2(y, x), Q.zero(y) & Q.zero(x)) == nan
+    assert refine(atan2(y, x), Q.zero(y) & Q.zero(x)) is nan
+
+
+def test_re():
+    assert refine(re(x), Q.real(x)) == x
+    assert refine(re(x), Q.imaginary(x)) == 0
+    assert refine(re(x+y), Q.real(x) & Q.real(y)) == x + y
+    assert refine(re(x+y), Q.real(x) & Q.imaginary(y)) == x
+    assert refine(re(x*y), Q.real(x) & Q.real(y)) == x * y
+    assert refine(re(x*y), Q.real(x) & Q.imaginary(y)) == 0
+    assert refine(re(x*y*z), Q.real(x) & Q.real(y) & Q.real(z)) == x * y * z
+
+
+def test_im():
+    assert refine(im(x), Q.imaginary(x)) == -I*x
+    assert refine(im(x), Q.real(x)) == 0
+    assert refine(im(x+y), Q.imaginary(x) & Q.imaginary(y)) == -I*x - I*y
+    assert refine(im(x+y), Q.real(x) & Q.imaginary(y)) == -I*y
+    assert refine(im(x*y), Q.imaginary(x) & Q.real(y)) == -I*x*y
+    assert refine(im(x*y), Q.imaginary(x) & Q.imaginary(y)) == 0
+    assert refine(im(1/x), Q.imaginary(x)) == -I/x
+    assert refine(im(x*y*z), Q.imaginary(x) & Q.imaginary(y)
+        & Q.imaginary(z)) == -I*x*y*z
+
+
+def test_complex():
+    assert refine(re(1/(x + I*y)), Q.real(x) & Q.real(y)) == \
+        x/(x**2 + y**2)
+    assert refine(im(1/(x + I*y)), Q.real(x) & Q.real(y)) == \
+        -y/(x**2 + y**2)
+    assert refine(re((w + I*x) * (y + I*z)), Q.real(w) & Q.real(x) & Q.real(y)
+        & Q.real(z)) == w*y - x*z
+    assert refine(im((w + I*x) * (y + I*z)), Q.real(w) & Q.real(x) & Q.real(y)
+        & Q.real(z)) == w*z + x*y
 
 
 def test_func_args():

sympy/assumptions/tests/test_satask.py

@@ -1,8 +1,8 @@
 from sympy.assumptions.satask import satask
 
-from sympy import symbols, Q, assuming, Implies, MatrixSymbol, I, pi, Rational
+from sympy import S, symbols, Q, assuming, Implies, MatrixSymbol, I, pi
 
-from sympy.utilities.pytest import raises, XFAIL, slow
+from sympy.utilities.pytest import raises, XFAIL
 
 
 x, y, z = symbols('x y z')
@@ -200,7 +200,7 @@ def test_odd_satask():
 
 def test_integer():
     assert satask(Q.integer(1)) is True
-    assert satask(Q.integer(Rational(1, 2))) is False
+    assert satask(Q.integer(S.Half)) is False
 
     assert satask(Q.integer(x + y), Q.integer(x) & Q.integer(y)) is True
     assert satask(Q.integer(x + y), Q.integer(x)) is None
@@ -264,7 +264,6 @@ def test_pos_neg():
     assert satask(Q.negative(x + y), Q.positive(x) & Q.positive(y)) is False
 
 
-@slow
 def test_pow_pos_neg():
     assert satask(Q.nonnegative(x**2), Q.positive(x)) is True
     assert satask(Q.nonpositive(x**2), Q.positive(x)) is False
@@ -323,3 +322,16 @@ def test_pow_pos_neg():
 
     # We could deduce things for negative powers if x is nonzero, but it
     # isn't implemented yet.
+
+
+def test_prime_composite():
+    assert satask(Q.prime(x), Q.composite(x)) is False
+    assert satask(Q.composite(x), Q.prime(x)) is False
+    assert satask(Q.composite(x), ~Q.prime(x)) is None
+    assert satask(Q.prime(x), ~Q.composite(x)) is None
+    # since 1 is neither prime nor composite the following should hold
+    assert satask(Q.prime(x), Q.integer(x) & Q.positive(x) & ~Q.composite(x)) is None
+    assert satask(Q.prime(2)) is True
+    assert satask(Q.prime(4)) is False
+    assert satask(Q.prime(1)) is False
+    assert satask(Q.composite(1)) is False

sympy/assumptions/tests/test_sathandlers.py

@@ -1,4 +1,5 @@
-from sympy import Mul, Basic, Q, Expr, And, symbols, Equivalent, Implies, Or
+from sympy import Mul, Basic, Q, Expr, And, symbols, Equivalent, Or
+from sympy.assumptions.cnf import to_NNF
 
 from sympy.assumptions.sathandlers import (ClassFactRegistry, AllArgs,
     UnevaluatedOnFree, AnyArgs, CheckOldAssump, ExactlyOneArg)
@@ -39,7 +40,7 @@ def test_UnevaluatedOnFree():
         Q.negative(y)))
 
     class MyUnevaluatedOnFree(UnevaluatedOnFree):
-        def apply(self):
+        def apply(self, expr=None):
             return self.args[0]
 
     a = MyUnevaluatedOnFree(Q.positive)
@@ -56,32 +57,32 @@ def test_UnevaluatedOnFree():
 def test_AllArgs():
     a = AllArgs(Q.zero)
     b = AllArgs(Q.positive | Q.negative)
-    assert a.rcall(x*y) == And(Q.zero(x), Q.zero(y))
-    assert b.rcall(x*y) == And(Q.positive(x) | Q.negative(x), Q.positive(y) | Q.negative(y))
+    assert a.rcall(x*y) == to_NNF(And(Q.zero(x), Q.zero(y)))
+    assert b.rcall(x*y) == to_NNF(And(Q.positive(x) | Q.negative(x), Q.positive(y) | Q.negative(y)))
 
 
 def test_AnyArgs():
     a = AnyArgs(Q.zero)
     b = AnyArgs(Q.positive & Q.negative)
-    assert a.rcall(x*y) == Or(Q.zero(x), Q.zero(y))
-    assert b.rcall(x*y) == Or(Q.positive(x) & Q.negative(x), Q.positive(y) & Q.negative(y))
+    assert a.rcall(x*y) == to_NNF(Or(Q.zero(x), Q.zero(y)))
+    assert b.rcall(x*y) == to_NNF(Or(Q.positive(x) & Q.negative(x), Q.positive(y) & Q.negative(y)))
 
 
 def test_CheckOldAssump():
     # TODO: Make these tests more complete
 
     class Test1(Expr):
-        def _eval_is_positive(self):
+        def _eval_is_extended_positive(self):
             return True
-        def _eval_is_negative(self):
+        def _eval_is_extended_negative(self):
             return False
 
     class Test2(Expr):
         def _eval_is_finite(self):
             return True
-        def _eval_is_positive(self):
+        def _eval_is_extended_positive(self):
             return True
-        def _eval_is_negative(self):
+        def _eval_is_extended_negative(self):
             return False
 
     t1 = Test1()
@@ -90,19 +91,19 @@ def test_CheckOldAssump():
     # We can't say if it's positive or negative in the old assumptions without
     # bounded. Remember, True means "no new knowledge", and
     # Q.positive(t2) means "t2 is positive."
-    assert CheckOldAssump(Q.positive(t1)) == True
-    assert CheckOldAssump(Q.negative(t1)) == ~Q.negative(t1)
+    assert CheckOldAssump(Q.positive(t1)) == to_NNF(True)
+    assert CheckOldAssump(Q.negative(t1)) == to_NNF(~Q.negative(t1))
 
-    assert CheckOldAssump(Q.positive(t2)) == Q.positive(t2)
-    assert CheckOldAssump(Q.negative(t2)) == ~Q.negative(t2)
+    assert CheckOldAssump(Q.positive(t2)) == to_NNF(Q.positive(t2))
+    assert CheckOldAssump(Q.negative(t2)) == to_NNF(~Q.negative(t2))
 
 
 def test_ExactlyOneArg():
     a = ExactlyOneArg(Q.zero)
     b = ExactlyOneArg(Q.positive | Q.negative)
-    assert a.rcall(x*y) == Or(Q.zero(x) & ~Q.zero(y), Q.zero(y) & ~Q.zero(x))
-    assert a.rcall(x*y*z) == Or(Q.zero(x) & ~Q.zero(y) & ~Q.zero(z), Q.zero(y)
-        & ~Q.zero(x) & ~Q.zero(z), Q.zero(z) & ~Q.zero(x) & ~Q.zero(y))
-    assert b.rcall(x*y) == Or((Q.positive(x) | Q.negative(x)) &
+    assert a.rcall(x*y) == to_NNF(Or(Q.zero(x) & ~Q.zero(y), Q.zero(y) & ~Q.zero(x)))
+    assert a.rcall(x*y*z) == to_NNF(Or(Q.zero(x) & ~Q.zero(y) & ~Q.zero(z), Q.zero(y)
+        & ~Q.zero(x) & ~Q.zero(z), Q.zero(z) & ~Q.zero(x) & ~Q.zero(y)))
+    assert b.rcall(x*y) == to_NNF(Or((Q.positive(x) | Q.negative(x)) &
         ~(Q.positive(y) | Q.negative(y)), (Q.positive(y) | Q.negative(y)) &
-        ~(Q.positive(x) | Q.negative(x)))
+        ~(Q.positive(x) | Q.negative(x))))

sympy/calculus/finite_diff.py

@@ -148,7 +148,7 @@ def finite_diff_weights(order, x_list, x0=S.One):
     derivative is wanted, weights for 0th, 1st and 2nd order are
     calculated "for free", so are formulae using subsets of ``x_list``.
     This is something one can take advantage of to save computational cost.
-    Be aware that one should define ``x_list`` from nearest to farest from
+    Be aware that one should define ``x_list`` from nearest to furthest from
     ``x0``. If not, subsets of ``x_list`` will yield poorer approximations,
     which might not grand an order of accuracy of ``len(x_list) - order``.
 
@@ -167,6 +167,9 @@ def finite_diff_weights(order, x_list, x0=S.One):
 
     """
     # The notation below closely corresponds to the one used in the paper.
+    order = S(order)
+    if not order.is_number:
+        raise ValueError("Cannot handle symbolic order.")
     if order < 0:
         raise ValueError("Negative derivative order illegal.")
     if int(order) != order:
@@ -175,10 +178,10 @@ def finite_diff_weights(order, x_list, x0=S.One):
     N = len(x_list) - 1
     delta = [[[0 for nu in range(N+1)] for n in range(N+1)] for
              m in range(M+1)]
-    delta[0][0][0] = S(1)
-    c1 = S(1)
+    delta[0][0][0] = S.One
+    c1 = S.One
     for n in range(1, N+1):
-        c2 = S(1)
+        c2 = S.One
         for nu in range(0, n):
             c3 = x_list[n]-x_list[nu]
             c2 = c2 * c3
@@ -195,7 +198,7 @@ def finite_diff_weights(order, x_list, x0=S.One):
     return delta
 
 
-def apply_finite_diff(order, x_list, y_list, x0=S(0)):
+def apply_finite_diff(order, x_list, y_list, x0=S.Zero):
     """
     Calculates the finite difference approximation of
     the derivative of requested order at ``x0`` from points
@@ -213,7 +216,7 @@ def apply_finite_diff(order, x_list, y_list, x0=S(0)):
         variable in x_list.
     x0: Number or Symbol
         At what value of the independent variable the derivative should be
-        evaluated. Defaults to S(0).
+        evaluated. Defaults to 0.
 
     Returns
     =======
@@ -387,6 +390,8 @@ def _as_finite_diff(derivative, points=1, x0=None, wrt=None):
         x0 = wrt
 
     if not iterable(points):
+        if getattr(points, 'is_Function', False) and wrt in points.args:
+            points = points.subs(wrt, x0)
         # points is simply the step-size, let's make it a
         # equidistant sequence centered around x0
         if order % 2 == 0:
@@ -413,6 +418,10 @@ as_finite_diff = deprecated(
     useinstead="Derivative.as_finite_difference",
     deprecated_since_version="1.1", issue=11410)(_as_finite_diff)
 
+as_finite_diff.__doc__ = """
+    Deprecated function. Use Diff.as_finite_difference instead.
+    """
+
 
 def differentiate_finite(expr, *symbols,
                          # points=1, x0=None, wrt=None, evaluate=True, #Py2:
@@ -465,7 +474,7 @@ def differentiate_finite(expr, *symbols,
     x0 = kwargs.pop('x0', None)
     wrt = kwargs.pop('wrt', None)
     evaluate = kwargs.pop('evaluate', False)
-    if kwargs != {}:
+    if kwargs:
         raise ValueError("Unknown kwargs: %s" % kwargs)
 
     Dexpr = expr.diff(*symbols, evaluate=evaluate)

sympy/calculus/tests/test_euler.py

@@ -11,6 +11,7 @@ def test_euler_interface():
     raises(TypeError, lambda: euler(D(x(t), t)*y(t), [x(t), y]))
     raises(ValueError, lambda: euler(D(x(t), t)*x(y), [x(t), x(y)]))
     raises(TypeError, lambda: euler(D(x(t), t)**2, x(0)))
+    raises(TypeError, lambda: euler(D(x(t), t)*y(t), [t]))
     assert euler(D(x(t), t)**2/2, {x(t)}) == [Eq(-D(x(t), t, t), 0)]
     assert euler(D(x(t), t)**2/2, x(t), {t}) == [Eq(-D(x(t), t, t), 0)]
 

sympy/calculus/tests/test_finite_diff.py

@@ -1,12 +1,12 @@
 from itertools import product
 
-from sympy import S, symbols, Function, exp
-from sympy.core.compatibility import range
-from sympy.utilities.pytest import warns_deprecated_sympy
+from sympy import S, symbols, Function, exp, diff, Rational
 from sympy.calculus.finite_diff import (
     apply_finite_diff, differentiate_finite, finite_diff_weights,
     as_finite_diff
 )
+from sympy.core.compatibility import range
+from sympy.utilities.pytest import raises, warns_deprecated_sympy
 
 
 def test_apply_finite_diff():
@@ -16,61 +16,62 @@ def test_apply_finite_diff():
             (f(x+h)-f(x-h))/(2*h)).simplify() == 0
 
     assert (apply_finite_diff(1, [5, 6, 7], [f(5), f(6), f(7)], 5) -
-            (-S(3)/2*f(5) + 2*f(6) - S(1)/2*f(7))).simplify() == 0
+            (Rational(-3, 2)*f(5) + 2*f(6) - S.Half*f(7))).simplify() == 0
+    raises(ValueError, lambda: apply_finite_diff(1, [x, h], [f(x)]))
 
 
 def test_finite_diff_weights():
 
     d = finite_diff_weights(1, [5, 6, 7], 5)
-    assert d[1][2] == [-S(3)/2, 2, -S(1)/2]
+    assert d[1][2] == [Rational(-3, 2), 2, Rational(-1, 2)]
 
     # Table 1, p. 702 in doi:10.1090/S0025-5718-1988-0935077-0
     # --------------------------------------------------------
     xl = [0, 1, -1, 2, -2, 3, -3, 4, -4]
 
     # d holds all coefficients
-    d = finite_diff_weights(4, xl, S(0))
+    d = finite_diff_weights(4, xl, S.Zero)
 
     # Zeroeth derivative
     for i in range(5):
-        assert d[0][i] == [S(1)] + [S(0)]*8
+        assert d[0][i] == [S.One] + [S.Zero]*8
 
     # First derivative
-    assert d[1][0] == [S(0)]*9
-    assert d[1][2] == [S(0), S(1)/2, -S(1)/2] + [S(0)]*6
-    assert d[1][4] == [S(0), S(2)/3, -S(2)/3, -S(1)/12, S(1)/12] + [S(0)]*4
-    assert d[1][6] == [S(0), S(3)/4, -S(3)/4, -S(3)/20, S(3)/20,
-                       S(1)/60, -S(1)/60] + [S(0)]*2
-    assert d[1][8] == [S(0), S(4)/5, -S(4)/5, -S(1)/5, S(1)/5,
-                       S(4)/105, -S(4)/105, -S(1)/280, S(1)/280]
+    assert d[1][0] == [S.Zero]*9
+    assert d[1][2] == [S.Zero, S.Half, Rational(-1, 2)] + [S.Zero]*6
+    assert d[1][4] == [S.Zero, Rational(2, 3), Rational(-2, 3), Rational(-1, 12), Rational(1, 12)] + [S.Zero]*4
+    assert d[1][6] == [S.Zero, Rational(3, 4), Rational(-3, 4), Rational(-3, 20), Rational(3, 20),
+                       Rational(1, 60), Rational(-1, 60)] + [S.Zero]*2
+    assert d[1][8] == [S.Zero, Rational(4, 5), Rational(-4, 5), Rational(-1, 5), Rational(1, 5),
+                       Rational(4, 105), Rational(-4, 105), Rational(-1, 280), Rational(1, 280)]
 
     # Second derivative
     for i in range(2):
-        assert d[2][i] == [S(0)]*9
-    assert d[2][2] == [-S(2), S(1), S(1)] + [S(0)]*6
-    assert d[2][4] == [-S(5)/2, S(4)/3, S(4)/3, -S(1)/12, -S(1)/12] + [S(0)]*4
-    assert d[2][6] == [-S(49)/18, S(3)/2, S(3)/2, -S(3)/20, -S(3)/20,
-                       S(1)/90, S(1)/90] + [S(0)]*2
-    assert d[2][8] == [-S(205)/72, S(8)/5, S(8)/5, -S(1)/5, -S(1)/5,
-                       S(8)/315, S(8)/315, -S(1)/560, -S(1)/560]
+        assert d[2][i] == [S.Zero]*9
+    assert d[2][2] == [-S(2), S.One, S.One] + [S.Zero]*6
+    assert d[2][4] == [Rational(-5, 2), Rational(4, 3), Rational(4, 3), Rational(-1, 12), Rational(-1, 12)] + [S.Zero]*4
+    assert d[2][6] == [Rational(-49, 18), Rational(3, 2), Rational(3, 2), Rational(-3, 20), Rational(-3, 20),
+                       Rational(1, 90), Rational(1, 90)] + [S.Zero]*2
+    assert d[2][8] == [Rational(-205, 72), Rational(8, 5), Rational(8, 5), Rational(-1, 5), Rational(-1, 5),
+                       Rational(8, 315), Rational(8, 315), Rational(-1, 560), Rational(-1, 560)]
 
     # Third derivative
     for i in range(3):
-        assert d[3][i] == [S(0)]*9
-    assert d[3][4] == [S(0), -S(1), S(1), S(1)/2, -S(1)/2] + [S(0)]*4
-    assert d[3][6] == [S(0), -S(13)/8, S(13)/8, S(1), -S(1),
-                       -S(1)/8, S(1)/8] + [S(0)]*2
-    assert d[3][8] == [S(0), -S(61)/30, S(61)/30, S(169)/120, -S(169)/120,
-                       -S(3)/10, S(3)/10, S(7)/240, -S(7)/240]
+        assert d[3][i] == [S.Zero]*9
+    assert d[3][4] == [S.Zero, -S.One, S.One, S.Half, Rational(-1, 2)] + [S.Zero]*4
+    assert d[3][6] == [S.Zero, Rational(-13, 8), Rational(13, 8), S.One, -S.One,
+                       Rational(-1, 8), Rational(1, 8)] + [S.Zero]*2
+    assert d[3][8] == [S.Zero, Rational(-61, 30), Rational(61, 30), Rational(169, 120), Rational(-169, 120),
+                       Rational(-3, 10), Rational(3, 10), Rational(7, 240), Rational(-7, 240)]
 
     # Fourth derivative
     for i in range(4):
-        assert d[4][i] == [S(0)]*9
-    assert d[4][4] == [S(6), -S(4), -S(4), S(1), S(1)] + [S(0)]*4
-    assert d[4][6] == [S(28)/3, -S(13)/2, -S(13)/2, S(2), S(2),
-                       -S(1)/6, -S(1)/6] + [S(0)]*2
-    assert d[4][8] == [S(91)/8, -S(122)/15, -S(122)/15, S(169)/60, S(169)/60,
-                       -S(2)/5, -S(2)/5, S(7)/240, S(7)/240]
+        assert d[4][i] == [S.Zero]*9
+    assert d[4][4] == [S(6), -S(4), -S(4), S.One, S.One] + [S.Zero]*4
+    assert d[4][6] == [Rational(28, 3), Rational(-13, 2), Rational(-13, 2), S(2), S(2),
+                       Rational(-1, 6), Rational(-1, 6)] + [S.Zero]*2
+    assert d[4][8] == [Rational(91, 8), Rational(-122, 15), Rational(-122, 15), Rational(169, 60), Rational(169, 60),
+                       Rational(-2, 5), Rational(-2, 5), Rational(7, 240), Rational(7, 240)]
 
     # Table 2, p. 703 in doi:10.1090/S0025-5718-1988-0935077-0
     # --------------------------------------------------------
@@ -82,47 +83,62 @@ def test_finite_diff_weights():
          i in range(4)]
 
     # Zeroth derivative
-    assert d[0][0][1] == [S(1)/2, S(1)/2]
-    assert d[1][0][3] == [-S(1)/16, S(9)/16, S(9)/16, -S(1)/16]
-    assert d[2][0][5] == [S(3)/256, -S(25)/256, S(75)/128, S(75)/128,
-                          -S(25)/256, S(3)/256]
-    assert d[3][0][7] == [-S(5)/2048, S(49)/2048, -S(245)/2048, S(1225)/2048,
-                          S(1225)/2048, -S(245)/2048, S(49)/2048, -S(5)/2048]
+    assert d[0][0][1] == [S.Half, S.Half]
+    assert d[1][0][3] == [Rational(-1, 16), Rational(9, 16), Rational(9, 16), Rational(-1, 16)]
+    assert d[2][0][5] == [Rational(3, 256), Rational(-25, 256), Rational(75, 128), Rational(75, 128),
+                          Rational(-25, 256), Rational(3, 256)]
+    assert d[3][0][7] == [Rational(-5, 2048), Rational(49, 2048), Rational(-245, 2048), Rational(1225, 2048),
+                          Rational(1225, 2048), Rational(-245, 2048), Rational(49, 2048), Rational(-5, 2048)]
 
     # First derivative
-    assert d[0][1][1] == [-S(1), S(1)]
-    assert d[1][1][3] == [S(1)/24, -S(9)/8, S(9)/8, -S(1)/24]
-    assert d[2][1][5] == [-S(3)/640, S(25)/384, -S(75)/64, S(75)/64,
-                          -S(25)/384, S(3)/640]
-    assert d[3][1][7] == [S(5)/7168, -S(49)/5120,  S(245)/3072, S(-1225)/1024,
-                          S(1225)/1024, -S(245)/3072, S(49)/5120, -S(5)/7168]
+    assert d[0][1][1] == [-S.One, S.One]
+    assert d[1][1][3] == [Rational(1, 24), Rational(-9, 8), Rational(9, 8), Rational(-1, 24)]
+    assert d[2][1][5] == [Rational(-3, 640), Rational(25, 384), Rational(-75, 64),
+                          Rational(75, 64), Rational(-25, 384), Rational(3, 640)]
+    assert d[3][1][7] == [Rational(5, 7168), Rational(-49, 5120),
+                          Rational(245, 3072), Rational(-1225, 1024),
+                          Rational(1225, 1024), Rational(-245, 3072),
+                          Rational(49, 5120), Rational(-5, 7168)]
 
     # Reasonably the rest of the table is also correct... (testing of that
     # deemed excessive at the moment)
+    raises(ValueError, lambda: finite_diff_weights(-1, [1, 2]))
+    raises(ValueError, lambda: finite_diff_weights(1.2, [1, 2]))
+    x = symbols('x')
+    raises(ValueError, lambda: finite_diff_weights(x, [1, 2]))
 
 
 def test_as_finite_diff():
     x = symbols('x')
     f = Function('f')
+    dx = Function('dx')
 
     with warns_deprecated_sympy():
         as_finite_diff(f(x).diff(x), [x-2, x-1, x, x+1, x+2])
 
+    # Use of undefined functions in ``points``
+    df_true = -f(x+dx(x)/2-dx(x+dx(x)/2)/2) / dx(x+dx(x)/2) \
+              + f(x+dx(x)/2+dx(x+dx(x)/2)/2) / dx(x+dx(x)/2)
+    df_test = diff(f(x), x).as_finite_difference(points=dx(x), x0=x+dx(x)/2)
+    assert (df_test - df_true).simplify() == 0
+
 
 def test_differentiate_finite():
     x, y = symbols('x y')
     f = Function('f')
     res0 = differentiate_finite(f(x, y) + exp(42), x, y, evaluate=True)
-    xm, xp, ym, yp = [v + sign*S(1)/2 for v, sign in product([x, y], [-1, 1])]
+    xm, xp, ym, yp = [v + sign*S.Half for v, sign in product([x, y], [-1, 1])]
     ref0 = f(xm, ym) + f(xp, yp) - f(xm, yp) - f(xp, ym)
     assert (res0 - ref0).simplify() == 0
 
     g = Function('g')
     res1 = differentiate_finite(f(x)*g(x) + 42, x, evaluate=True)
-    ref1 = (-f(x - S(1)/2) + f(x + S(1)/2))*g(x) + \
-           (-g(x - S(1)/2) + g(x + S(1)/2))*f(x)
+    ref1 = (-f(x - S.Half) + f(x + S.Half))*g(x) + \
+           (-g(x - S.Half) + g(x + S.Half))*f(x)
     assert (res1 - ref1).simplify() == 0
 
     res2 = differentiate_finite(f(x) + x**3 + 42, x, points=[x-1, x+1])
     ref2 = (f(x + 1) + (x + 1)**3 - f(x - 1) - (x - 1)**3)/2
     assert (res2 - ref2).simplify() == 0
+    raises(ValueError, lambda: differentiate_finite(f(x)*g(x), x,
+                                                    pints=[x-1, x+1]))

sympy/calculus/tests/test_singularities.py

@@ -1,4 +1,4 @@
-from sympy import Symbol, exp, log, oo, S, I, sqrt
+from sympy import Symbol, exp, log, oo, S, I, sqrt, Rational
 from sympy.calculus.singularities import (
     singularities,
     is_increasing,
@@ -8,7 +8,7 @@ from sympy.calculus.singularities import (
     is_monotonic
 )
 from sympy.sets import Interval, FiniteSet
-from sympy.utilities.pytest import XFAIL
+from sympy.utilities.pytest import XFAIL, raises
 from sympy.abc import x, y
 
 
@@ -65,7 +65,7 @@ def test_is_decreasing():
 
     assert is_decreasing(1/(x**2 - 3*x), Interval.open(1.5, 3))
     assert is_decreasing(1/(x**2 - 3*x), Interval.Lopen(3, oo))
-    assert not is_decreasing(1/(x**2 - 3*x), Interval.Ropen(-oo, S(3)/2))
+    assert not is_decreasing(1/(x**2 - 3*x), Interval.Ropen(-oo, Rational(3, 2)))
     assert not is_decreasing(-x**2, Interval(-oo, 0))
     assert not is_decreasing(-x**2*b, Interval(-oo, 0), x)
 
@@ -74,7 +74,7 @@ def test_is_strictly_decreasing():
     """Test whether is_strictly_decreasing returns correct value."""
     assert is_strictly_decreasing(1/(x**2 - 3*x), Interval.Lopen(3, oo))
     assert not is_strictly_decreasing(
-        1/(x**2 - 3*x), Interval.Ropen(-oo, S(3)/2))
+        1/(x**2 - 3*x), Interval.Ropen(-oo, Rational(3, 2)))
     assert not is_strictly_decreasing(-x**2, Interval(-oo, 0))
     assert not is_strictly_decreasing(1)
     assert is_strictly_decreasing(1/(x**2 - 3*x), Interval.open(1.5, 3))
@@ -87,3 +87,4 @@ def test_is_monotonic():
     assert is_monotonic(x**3 - 3*x**2 + 4*x, S.Reals)
     assert not is_monotonic(-x**2, S.Reals)
     assert is_monotonic(x**2 + y + 1, Interval(1, 2), x)
+    raises(NotImplementedError, lambda: is_monotonic(x**2 + y + 1))

sympy/calculus/tests/test_util.py

@@ -1,12 +1,12 @@
 from sympy import (Symbol, S, exp, log, sqrt, oo, E, zoo, pi, tan, sin, cos,
-                   cot, sec, csc, Abs, symbols, I, re, Lambda, simplify,
-                   ImageSet)
+                   cot, sec, csc, Abs, symbols, I, re, simplify,
+                   expint, Rational)
 from sympy.calculus.util import (function_range, continuous_domain, not_empty_in,
                                  periodicity, lcim, AccumBounds, is_convex,
                                  stationary_points, minimum, maximum)
 from sympy.core import Add, Mul, Pow
 from sympy.sets.sets import (Interval, FiniteSet, EmptySet, Complement,
-                            Union, Intersection)
+                            Union)
 from sympy.utilities.pytest import raises
 from sympy.abc import x
 
@@ -24,7 +24,7 @@ def test_function_range():
     assert function_range(tan(x), x, Interval(pi/2, pi)
         ) == Interval(-oo, 0)
     assert function_range((x + 3)/(x - 2), x, Interval(-5, 5)
-        ) == Union(Interval(-oo, S(2)/7), Interval(S(8)/3, oo))
+        ) == Union(Interval(-oo, Rational(2, 7)), Interval(Rational(8, 3), oo))
     assert function_range(1/(x**2), x, Interval(-1, 1)
         ) == Interval(1, oo)
     assert function_range(exp(x), x, Interval(-1, 1)
@@ -41,8 +41,11 @@ def test_function_range():
         ) == Union(Interval(-sin(3), 1), FiniteSet(sin(4)))
     assert function_range(cos(x), x, Interval(-oo, -4)
         ) == Interval(-1, 1)
+    assert function_range(cos(x), x, S.EmptySet) == S.EmptySet
     raises(NotImplementedError, lambda : function_range(
         exp(x)*(sin(x) - cos(x))/2 - x, x, S.Reals))
+    raises(NotImplementedError, lambda : function_range(
+        sin(x) + x, x, S.Reals)) # issue 13273
     raises(NotImplementedError, lambda : function_range(
         log(x), x, S.Integers))
     raises(NotImplementedError, lambda : function_range(
@@ -53,12 +56,12 @@ def test_continuous_domain():
     x = Symbol('x')
     assert continuous_domain(sin(x), x, Interval(0, 2*pi)) == Interval(0, 2*pi)
     assert continuous_domain(tan(x), x, Interval(0, 2*pi)) == \
-        Union(Interval(0, pi/2, False, True), Interval(pi/2, 3*pi/2, True, True),
-              Interval(3*pi/2, 2*pi, True, False))
+        Union(Interval(0, pi/2, False, True), Interval(pi/2, pi*Rational(3, 2), True, True),
+              Interval(pi*Rational(3, 2), 2*pi, True, False))
     assert continuous_domain((x - 1)/((x - 1)**2), x, S.Reals) == \
         Union(Interval(-oo, 1, True, True), Interval(1, oo, True, True))
     assert continuous_domain(log(x) + log(4*x - 1), x, S.Reals) == \
-        Interval(S(1)/4, oo, True, True)
+        Interval(Rational(1, 4), oo, True, True)
     assert continuous_domain(1/sqrt(x - 3), x, S.Reals) == Interval(3, oo, True, True)
     assert continuous_domain(1/x - 2, x, S.Reals) == \
         Union(Interval.open(-oo, 0), Interval.open(0, oo))
@@ -68,12 +71,12 @@ def test_continuous_domain():
 
 def test_not_empty_in():
     assert not_empty_in(FiniteSet(x, 2*x).intersect(Interval(1, 2, True, False)), x) == \
-        Interval(S(1)/2, 2, True, False)
+        Interval(S.Half, 2, True, False)
     assert not_empty_in(FiniteSet(x, x**2).intersect(Interval(1, 2)), x) == \
         Union(Interval(-sqrt(2), -1), Interval(1, 2))
     assert not_empty_in(FiniteSet(x**2 + x, x).intersect(Interval(2, 4)), x) == \
-        Union(Interval(-sqrt(17)/2 - S(1)/2, -2),
-              Interval(1, -S(1)/2 + sqrt(17)/2), Interval(2, 4))
+        Union(Interval(-sqrt(17)/2 - S.Half, -2),
+              Interval(1, Rational(-1, 2) + sqrt(17)/2), Interval(2, 4))
     assert not_empty_in(FiniteSet(x/(x - 1)).intersect(S.Reals), x) == \
         Complement(S.Reals, FiniteSet(1))
     assert not_empty_in(FiniteSet(a/(a - 1)).intersect(S.Reals), a) == \
@@ -81,7 +84,9 @@ def test_not_empty_in():
     assert not_empty_in(FiniteSet((x**2 - 3*x + 2)/(x - 1)).intersect(S.Reals), x) == \
         Complement(S.Reals, FiniteSet(1))
     assert not_empty_in(FiniteSet(3, 4, x/(x - 1)).intersect(Interval(2, 3)), x) == \
-        Union(Interval(S(3)/2, 2), FiniteSet(3))
+        Interval(-oo, oo)
+    assert not_empty_in(FiniteSet(4, x/(x - 1)).intersect(Interval(2, 3)), x) == \
+        Interval(S(3)/2, 2)
     assert not_empty_in(FiniteSet(x/(x**2 - 1)).intersect(S.Reals), x) == \
         Complement(S.Reals, FiniteSet(-1, 1))
     assert not_empty_in(FiniteSet(x, x**2).intersect(Union(Interval(1, 3, True, True),
@@ -91,6 +96,10 @@ def test_not_empty_in():
     assert not_empty_in(FiniteSet(1).intersect(Interval(3, 4)), x) == S.EmptySet
     assert not_empty_in(FiniteSet(x**2/(x + 2)).intersect(Interval(1, oo)), x) == \
         Union(Interval(-2, -1, True, False), Interval(2, oo))
+    raises(ValueError, lambda: not_empty_in(x))
+    raises(ValueError, lambda: not_empty_in(Interval(0, 1), x))
+    raises(NotImplementedError,
+           lambda: not_empty_in(FiniteSet(x).intersect(S.Reals), x, a))
 
 
 def test_periodicity():
@@ -118,12 +127,12 @@ def test_periodicity():
     assert periodicity(2*tan(x)**2, x) == pi
     assert periodicity(sin(x%4), x) == 4
     assert periodicity(sin(x)%4, x) == 2*pi
-    assert periodicity(tan((3*x-2)%4), x) == S(4)/3
+    assert periodicity(tan((3*x-2)%4), x) == Rational(4, 3)
     assert periodicity((sqrt(2)*(x+1)+x) % 3, x) == 3 / (sqrt(2)+1)
-    assert periodicity((x**2+1) % x, x) == None
+    assert periodicity((x**2+1) % x, x) is None
     assert periodicity(sin(re(x)), x) == 2*pi
-    assert periodicity(sin(x)**2 + cos(x)**2, x) == S.Zero
-    assert periodicity(tan(x), y) == S.Zero
+    assert periodicity(sin(x)**2 + cos(x)**2, x) is S.Zero
+    assert periodicity(tan(x), y) is S.Zero
     assert periodicity(sin(x) + I*cos(x), x) == 2*pi
     assert periodicity(x - sin(2*y), y) == pi
 
@@ -156,6 +165,9 @@ def test_periodicity():
     assert periodicity((x**2 + 4)%2, x) is None
     assert periodicity((E**x)%3, x) is None
 
+    assert periodicity(sin(expint(1, x))/expint(1, x), x) is None
+
+
 def test_periodicity_check():
     x = Symbol('x')
     y = Symbol('y')
@@ -170,18 +182,20 @@ def test_periodicity_check():
 def test_lcim():
     from sympy import pi
 
-    assert lcim([S(1)/2, S(2), S(3)]) == 6
+    assert lcim([S.Half, S(2), S(3)]) == 6
     assert lcim([pi/2, pi/4, pi]) == pi
     assert lcim([2*pi, pi/2]) == 2*pi
-    assert lcim([S(1), 2*pi]) is None
-    assert lcim([S(2) + 2*E, E/3 + S(1)/3, S(1) + E]) == S(2) + 2*E
+    assert lcim([S.One, 2*pi]) is None
+    assert lcim([S(2) + 2*E, E/3 + Rational(1, 3), S.One + E]) == S(2) + 2*E
+
 
 def test_is_convex():
-
     assert is_convex(1/x, x, domain=Interval(0, oo)) == True
     assert is_convex(1/x, x, domain=Interval(-oo, 0)) == False
     assert is_convex(x**2, x, domain=Interval(0, oo)) == True
     assert is_convex(log(x), x) == False
+    raises(NotImplementedError, lambda: is_convex(log(x), x, a))
+
 
 def test_stationary_points():
     x, y = symbols('x y')
@@ -193,11 +207,11 @@ def test_stationary_points():
     assert stationary_points(tan(x), x,
         ) == EmptySet()
     assert stationary_points(sin(x)*cos(x), x, Interval(0, pi)
-        ) == {pi/4, 3*pi/4}
+        ) == {pi/4, pi*Rational(3, 4)}
     assert stationary_points(sec(x), x, Interval(0, pi)
         ) == {0, pi}
     assert stationary_points((x+3)*(x-2), x
-        ) == FiniteSet(-S.Half)
+        ) == FiniteSet(Rational(-1, 2))
     assert stationary_points((x + 3)/(x - 2), x, Interval(-5, 5)
         ) == EmptySet()
     assert stationary_points((x**2+3)/(x-2), x
@@ -205,7 +219,7 @@ def test_stationary_points():
     assert stationary_points((x**2+3)/(x-2), x, Interval(0, 5)
         ) == {2 + sqrt(7)}
     assert stationary_points(x**4 + x**3 - 5*x**2, x, S.Reals
-        ) == FiniteSet(-2, 0, S(5)/4)
+        ) == FiniteSet(-2, 0, Rational(5, 4))
     assert stationary_points(exp(x), x
         ) == EmptySet()
     assert stationary_points(log(x) - x, x, S.Reals
@@ -214,25 +228,27 @@ def test_stationary_points():
         ) == {0, -pi, pi}
     assert stationary_points(y, x, S.Reals
         ) == S.Reals
+    assert stationary_points(y, x, S.EmptySet) == S.EmptySet
+
 
 def test_maximum():
     x, y = symbols('x y')
-    assert maximum(sin(x), x) == S.One
+    assert maximum(sin(x), x) is S.One
     assert maximum(sin(x), x, Interval(0, 1)) == sin(1)
-    assert maximum(tan(x), x) == oo
-    assert maximum(tan(x), x, Interval(-pi/4, pi/4)) == S.One
+    assert maximum(tan(x), x) is oo
+    assert maximum(tan(x), x, Interval(-pi/4, pi/4)) is S.One
     assert maximum(sin(x)*cos(x), x, S.Reals) == S.Half
-    assert simplify(maximum(sin(x)*cos(x), x, Interval(3*pi/8, 5*pi/8))
+    assert simplify(maximum(sin(x)*cos(x), x, Interval(pi*Rational(3, 8), pi*Rational(5, 8)))
         ) == sqrt(2)/4
-    assert maximum((x+3)*(x-2), x) == oo
+    assert maximum((x+3)*(x-2), x) is oo
     assert maximum((x+3)*(x-2), x, Interval(-5, 0)) == S(14)
-    assert maximum((x+3)/(x-2), x, Interval(-5, 0)) == S(2)/7
+    assert maximum((x+3)/(x-2), x, Interval(-5, 0)) == Rational(2, 7)
     assert simplify(maximum(-x**4-x**3+x**2+10, x)
-        ) == 41*sqrt(41)/512 + S(5419)/512
+        ) == 41*sqrt(41)/512 + Rational(5419, 512)
     assert maximum(exp(x), x, Interval(-oo, 2)) == exp(2)
-    assert maximum(log(x) - x, x, S.Reals) == -S.One
+    assert maximum(log(x) - x, x, S.Reals) is S.NegativeOne
     assert maximum(cos(x), x, Union(Interval(0, 5), Interval(-6, -3))
-        ) == S.One
+        ) is S.One
     assert maximum(cos(x)-sin(x), x, S.Reals) == sqrt(2)
     assert maximum(y, x, S.Reals) == y
 
@@ -241,25 +257,26 @@ def test_maximum():
     raises(ValueError, lambda : maximum(1/(x**2 + y**2 + 1), x, S.EmptySet))
     raises(ValueError, lambda : maximum(sin(x), sin(x)))
     raises(ValueError, lambda : maximum(sin(x), x*y, S.EmptySet))
-    raises(ValueError, lambda : maximum(sin(x), S(1)))
+    raises(ValueError, lambda : maximum(sin(x), S.One))
+
 
 def test_minimum():
     x, y = symbols('x y')
 
-    assert minimum(sin(x), x) == -S.One
+    assert minimum(sin(x), x) is S.NegativeOne
     assert minimum(sin(x), x, Interval(1, 4)) == sin(4)
-    assert minimum(tan(x), x) == -oo
-    assert minimum(tan(x), x, Interval(-pi/4, pi/4)) == -S.One
-    assert minimum(sin(x)*cos(x), x, S.Reals) == -S.Half
-    assert simplify(minimum(sin(x)*cos(x), x, Interval(3*pi/8, 5*pi/8))
+    assert minimum(tan(x), x) is -oo
+    assert minimum(tan(x), x, Interval(-pi/4, pi/4)) is S.NegativeOne
+    assert minimum(sin(x)*cos(x), x, S.Reals) == Rational(-1, 2)
+    assert simplify(minimum(sin(x)*cos(x), x, Interval(pi*Rational(3, 8), pi*Rational(5, 8)))
         ) == -sqrt(2)/4
-    assert minimum((x+3)*(x-2), x) == -S(25)/4
-    assert minimum((x+3)/(x-2), x, Interval(-5, 0)) == -S(3)/2
+    assert minimum((x+3)*(x-2), x) == Rational(-25, 4)
+    assert minimum((x+3)/(x-2), x, Interval(-5, 0)) == Rational(-3, 2)
     assert minimum(x**4-x**3+x**2+10, x) == S(10)
     assert minimum(exp(x), x, Interval(-2, oo)) == exp(-2)
-    assert minimum(log(x) - x, x, S.Reals) == -oo
+    assert minimum(log(x) - x, x, S.Reals) is -oo
     assert minimum(cos(x), x, Union(Interval(0, 5), Interval(-6, -3))
-        ) == -S.One
+        ) is S.NegativeOne
     assert minimum(cos(x)-sin(x), x, S.Reals) == -sqrt(2)
     assert minimum(y, x, S.Reals) == y
 
@@ -268,15 +285,16 @@ def test_minimum():
     raises(ValueError, lambda : minimum(1/(x**2 + y**2 + 1), x, S.EmptySet))
     raises(ValueError, lambda : minimum(sin(x), sin(x)))
     raises(ValueError, lambda : minimum(sin(x), x*y, S.EmptySet))
-    raises(ValueError, lambda : minimum(sin(x), S(1)))
+    raises(ValueError, lambda : minimum(sin(x), S.One))
+
 
 def test_AccumBounds():
     assert AccumBounds(1, 2).args == (1, 2)
-    assert AccumBounds(1, 2).delta == S(1)
-    assert AccumBounds(1, 2).mid == S(3)/2
+    assert AccumBounds(1, 2).delta is S.One
+    assert AccumBounds(1, 2).mid == Rational(3, 2)
     assert AccumBounds(1, 3).is_real == True
 
-    assert AccumBounds(1, 1) == S(1)
+    assert AccumBounds(1, 1) is S.One
 
     assert AccumBounds(1, 2) + 1 == AccumBounds(2, 3)
     assert 1 + AccumBounds(1, 2) == AccumBounds(2, 3)
@@ -293,24 +311,27 @@ def test_AccumBounds():
     assert AccumBounds(1, 2) - x == Add(AccumBounds(1, 2), -x)
 
     assert AccumBounds(-oo, 1) + oo == AccumBounds(-oo, oo)
-    assert AccumBounds(1, oo) + oo == oo
+    assert AccumBounds(1, oo) + oo is oo
     assert AccumBounds(1, oo) - oo == AccumBounds(-oo, oo)
-    assert (-oo - AccumBounds(-1, oo)) == -oo
-    assert AccumBounds(-oo, 1) - oo == -oo
+    assert (-oo - AccumBounds(-1, oo)) is -oo
+    assert AccumBounds(-oo, 1) - oo is -oo
 
     assert AccumBounds(1, oo) - oo == AccumBounds(-oo, oo)
     assert AccumBounds(-oo, 1) - (-oo) == AccumBounds(-oo, oo)
     assert (oo - AccumBounds(1, oo)) == AccumBounds(-oo, oo)
-    assert (-oo - AccumBounds(1, oo)) == -oo
+    assert (-oo - AccumBounds(1, oo)) is -oo
 
-    assert AccumBounds(1, 2)/2 == AccumBounds(S(1)/2, 1)
-    assert 2/AccumBounds(2, 3) == AccumBounds(S(2)/3, 1)
+    assert AccumBounds(1, 2)/2 == AccumBounds(S.Half, 1)
+    assert 2/AccumBounds(2, 3) == AccumBounds(Rational(2, 3), 1)
     assert 1/AccumBounds(-1, 1) == AccumBounds(-oo, oo)
 
     assert abs(AccumBounds(1, 2)) == AccumBounds(1, 2)
     assert abs(AccumBounds(-2, -1)) == AccumBounds(1, 2)
     assert abs(AccumBounds(-2, 1)) == AccumBounds(0, 2)
     assert abs(AccumBounds(-1, 2)) == AccumBounds(0, 2)
+    c = Symbol('c')
+    raises(ValueError, lambda: AccumBounds(0, c))
+    raises(ValueError, lambda: AccumBounds(1, -1))
 
 
 def test_AccumBounds_mul():
@@ -335,26 +356,26 @@ def test_AccumBounds_mul():
 
 
 def test_AccumBounds_div():
-    assert AccumBounds(-1, 3)/AccumBounds(3, 4) == AccumBounds(-S(1)/3, 1)
+    assert AccumBounds(-1, 3)/AccumBounds(3, 4) == AccumBounds(Rational(-1, 3), 1)
     assert AccumBounds(-2, 4)/AccumBounds(-3, 4) == AccumBounds(-oo, oo)
-    assert AccumBounds(-3, -2)/AccumBounds(-4, 0) == AccumBounds(S(1)/2, oo)
+    assert AccumBounds(-3, -2)/AccumBounds(-4, 0) == AccumBounds(S.Half, oo)
 
     # these two tests can have a better answer
     # after Union of AccumBounds is improved
     assert AccumBounds(-3, -2)/AccumBounds(-2, 1) == AccumBounds(-oo, oo)
     assert AccumBounds(2, 3)/AccumBounds(-2, 2) == AccumBounds(-oo, oo)
 
-    assert AccumBounds(-3, -2)/AccumBounds(0, 4) == AccumBounds(-oo, -S(1)/2)
-    assert AccumBounds(2, 4)/AccumBounds(-3, 0) == AccumBounds(-oo, -S(2)/3)
-    assert AccumBounds(2, 4)/AccumBounds(0, 3) == AccumBounds(S(2)/3, oo)
+    assert AccumBounds(-3, -2)/AccumBounds(0, 4) == AccumBounds(-oo, Rational(-1, 2))
+    assert AccumBounds(2, 4)/AccumBounds(-3, 0) == AccumBounds(-oo, Rational(-2, 3))
+    assert AccumBounds(2, 4)/AccumBounds(0, 3) == AccumBounds(Rational(2, 3), oo)
 
     assert AccumBounds(0, 1)/AccumBounds(0, 1) == AccumBounds(0, oo)
     assert AccumBounds(-1, 0)/AccumBounds(0, 1) == AccumBounds(-oo, 0)
     assert AccumBounds(-1, 2)/AccumBounds(-2, 2) == AccumBounds(-oo, oo)
 
     assert 1/AccumBounds(-1, 2) == AccumBounds(-oo, oo)
-    assert 1/AccumBounds(0, 2) == AccumBounds(S(1)/2, oo)
-    assert (-1)/AccumBounds(0, 2) == AccumBounds(-oo, -S(1)/2)
+    assert 1/AccumBounds(0, 2) == AccumBounds(S.Half, oo)
+    assert (-1)/AccumBounds(0, 2) == AccumBounds(-oo, Rational(-1, 2))
     assert 1/AccumBounds(-oo, 0) == AccumBounds(-oo, 0)
     assert 1/AccumBounds(-1, 0) == AccumBounds(-oo, -1)
     assert (-2)/AccumBounds(-oo, 0) == AccumBounds(0, oo)
@@ -389,51 +410,51 @@ def test_AccumBounds_pow():
     assert AccumBounds(-1, 2)**3 == AccumBounds(-1, 8)
     assert AccumBounds(-1, 1)**0 == 1
 
-    assert AccumBounds(1, 2)**(S(5)/2) == AccumBounds(1, 4*sqrt(2))
-    assert AccumBounds(-1, 2)**(S(1)/3) == AccumBounds(-1, 2**(S(1)/3))
-    assert AccumBounds(0, 2)**(S(1)/2) == AccumBounds(0, sqrt(2))
+    assert AccumBounds(1, 2)**Rational(5, 2) == AccumBounds(1, 4*sqrt(2))
+    assert AccumBounds(-1, 2)**Rational(1, 3) == AccumBounds(-1, 2**Rational(1, 3))
+    assert AccumBounds(0, 2)**S.Half == AccumBounds(0, sqrt(2))
 
-    assert AccumBounds(-4, 2)**(S(2)/3) == AccumBounds(0, 2*2**(S(1)/3))
+    assert AccumBounds(-4, 2)**Rational(2, 3) == AccumBounds(0, 2*2**Rational(1, 3))
 
-    assert AccumBounds(-1, 5)**(S(1)/2) == AccumBounds(0, sqrt(5))
-    assert AccumBounds(-oo, 2)**(S(1)/2) == AccumBounds(0, sqrt(2))
-    assert AccumBounds(-2, 3)**(S(-1)/4) == AccumBounds(0, oo)
+    assert AccumBounds(-1, 5)**S.Half == AccumBounds(0, sqrt(5))
+    assert AccumBounds(-oo, 2)**S.Half == AccumBounds(0, sqrt(2))
+    assert AccumBounds(-2, 3)**Rational(-1, 4) == AccumBounds(0, oo)
 
-    assert AccumBounds(1, 5)**(-2) == AccumBounds(S(1)/25, 1)
+    assert AccumBounds(1, 5)**(-2) == AccumBounds(Rational(1, 25), 1)
     assert AccumBounds(-1, 3)**(-2) == AccumBounds(0, oo)
-    assert AccumBounds(0, 2)**(-2) == AccumBounds(S(1)/4, oo)
+    assert AccumBounds(0, 2)**(-2) == AccumBounds(Rational(1, 4), oo)
     assert AccumBounds(-1, 2)**(-3) == AccumBounds(-oo, oo)
-    assert AccumBounds(-3, -2)**(-3) == AccumBounds(S(-1)/8, -S(1)/27)
-    assert AccumBounds(-3, -2)**(-2) == AccumBounds(S(1)/9, S(1)/4)
-    assert AccumBounds(0, oo)**(S(1)/2) == AccumBounds(0, oo)
-    assert AccumBounds(-oo, -1)**(S(1)/3) == AccumBounds(-oo, -1)
-    assert AccumBounds(-2, 3)**(-S(1)/3) == AccumBounds(-oo, oo)
+    assert AccumBounds(-3, -2)**(-3) == AccumBounds(Rational(-1, 8), Rational(-1, 27))
+    assert AccumBounds(-3, -2)**(-2) == AccumBounds(Rational(1, 9), Rational(1, 4))
+    assert AccumBounds(0, oo)**S.Half == AccumBounds(0, oo)
+    assert AccumBounds(-oo, -1)**Rational(1, 3) == AccumBounds(-oo, -1)
+    assert AccumBounds(-2, 3)**(Rational(-1, 3)) == AccumBounds(-oo, oo)
     assert AccumBounds(-oo, 0)**(-2) == AccumBounds(0, oo)
-    assert AccumBounds(-2, 0)**(-2) == AccumBounds(S(1)/4, oo)
+    assert AccumBounds(-2, 0)**(-2) == AccumBounds(Rational(1, 4), oo)
 
-    assert AccumBounds(S(1)/3, S(1)/2)**oo == S(0)
-    assert AccumBounds(0, S(1)/2)**oo == S(0)
-    assert AccumBounds(S(1)/2, 1)**oo == AccumBounds(0, oo)
+    assert AccumBounds(Rational(1, 3), S.Half)**oo is S.Zero
+    assert AccumBounds(0, S.Half)**oo is S.Zero
+    assert AccumBounds(S.Half, 1)**oo == AccumBounds(0, oo)
     assert AccumBounds(0, 1)**oo == AccumBounds(0, oo)
-    assert AccumBounds(2, 3)**oo == oo
+    assert AccumBounds(2, 3)**oo is oo
     assert AccumBounds(1, 2)**oo == AccumBounds(0, oo)
-    assert AccumBounds(S(1)/2, 3)**oo == AccumBounds(0, oo)
-    assert AccumBounds(-S(1)/3, -S(1)/4)**oo == S(0)
-    assert AccumBounds(-1, -S(1)/2)**oo == AccumBounds(-oo, oo)
+    assert AccumBounds(S.Half, 3)**oo == AccumBounds(0, oo)
+    assert AccumBounds(Rational(-1, 3), Rational(-1, 4))**oo is S.Zero
+    assert AccumBounds(-1, Rational(-1, 2))**oo == AccumBounds(-oo, oo)
     assert AccumBounds(-3, -2)**oo == FiniteSet(-oo, oo)
     assert AccumBounds(-2, -1)**oo == AccumBounds(-oo, oo)
-    assert AccumBounds(-2, -S(1)/2)**oo == AccumBounds(-oo, oo)
-    assert AccumBounds(-S(1)/2, S(1)/2)**oo == S(0)
-    assert AccumBounds(-S(1)/2, 1)**oo == AccumBounds(0, oo)
-    assert AccumBounds(-S(2)/3, 2)**oo == AccumBounds(0, oo)
+    assert AccumBounds(-2, Rational(-1, 2))**oo == AccumBounds(-oo, oo)
+    assert AccumBounds(Rational(-1, 2), S.Half)**oo is S.Zero
+    assert AccumBounds(Rational(-1, 2), 1)**oo == AccumBounds(0, oo)
+    assert AccumBounds(Rational(-2, 3), 2)**oo == AccumBounds(0, oo)
     assert AccumBounds(-1, 1)**oo == AccumBounds(-oo, oo)
-    assert AccumBounds(-1, S(1)/2)**oo == AccumBounds(-oo, oo)
+    assert AccumBounds(-1, S.Half)**oo == AccumBounds(-oo, oo)
     assert AccumBounds(-1, 2)**oo == AccumBounds(-oo, oo)
-    assert AccumBounds(-2, S(1)/2)**oo == AccumBounds(-oo, oo)
+    assert AccumBounds(-2, S.Half)**oo == AccumBounds(-oo, oo)
 
     assert AccumBounds(1, 2)**x == Pow(AccumBounds(1, 2), x, evaluate=False)
 
-    assert AccumBounds(2, 3)**(-oo) == S(0)
+    assert AccumBounds(2, 3)**(-oo) is S.Zero
     assert AccumBounds(0, 2)**(-oo) == AccumBounds(0, oo)
     assert AccumBounds(-1, 2)**(-oo) == AccumBounds(-oo, oo)
 
@@ -456,9 +477,25 @@ def test_comparison_AccumBounds():
     assert (AccumBounds(1, 3) < AccumBounds(2, 4)).rel_op == '<'
     assert (AccumBounds(1, 3) < AccumBounds(-2, 0)) == S.false
 
+    assert (AccumBounds(1, 3) <= AccumBounds(4, 6)) == S.true
+    assert (AccumBounds(1, 3) <= AccumBounds(-2, 0)) == S.false
+
+    assert (AccumBounds(1, 3) > AccumBounds(4, 6)) == S.false
+    assert (AccumBounds(1, 3) > AccumBounds(-2, 0)) == S.true
+
+    assert (AccumBounds(1, 3) >= AccumBounds(4, 6)) == S.false
+    assert (AccumBounds(1, 3) >= AccumBounds(-2, 0)) == S.true
+
     # issue 13499
     assert (cos(x) > 0).subs(x, oo) == (AccumBounds(-1, 1) > 0)
 
+    c = Symbol('c')
+    raises(TypeError, lambda: (AccumBounds(0, 1) < c))
+    raises(TypeError, lambda: (AccumBounds(0, 1) <= c))
+    raises(TypeError, lambda: (AccumBounds(0, 1) > c))
+    raises(TypeError, lambda: (AccumBounds(0, 1) >= c))
+
+
 def test_contains_AccumBounds():
     assert (1 in AccumBounds(1, 2)) == S.true
     raises(TypeError, lambda: a in AccumBounds(1, 2))
@@ -473,3 +510,26 @@ def test_contains_AccumBounds():
     import itertools
     for perm in itertools.permutations([0, AccumBounds(-1, 1), x]):
         assert Mul(*perm) == 0
+
+
+def test_intersection_AccumBounds():
+    assert AccumBounds(0, 3).intersection(AccumBounds(1, 2)) == AccumBounds(1, 2)
+    assert AccumBounds(0, 3).intersection(AccumBounds(1, 4)) == AccumBounds(1, 3)
+    assert AccumBounds(0, 3).intersection(AccumBounds(-1, 2)) == AccumBounds(0, 2)
+    assert AccumBounds(0, 3).intersection(AccumBounds(-1, 4)) == AccumBounds(0, 3)
+    assert AccumBounds(0, 1).intersection(AccumBounds(2, 3)) == S.EmptySet
+    raises(TypeError, lambda: AccumBounds(0, 3).intersection(1))
+
+
+def test_union_AccumBounds():
+    assert AccumBounds(0, 3).union(AccumBounds(1, 2)) == AccumBounds(0, 3)
+    assert AccumBounds(0, 3).union(AccumBounds(1, 4)) == AccumBounds(0, 4)
+    assert AccumBounds(0, 3).union(AccumBounds(-1, 2)) == AccumBounds(-1, 3)
+    assert AccumBounds(0, 3).union(AccumBounds(-1, 4)) == AccumBounds(-1, 4)
+    raises(TypeError, lambda: AccumBounds(0, 3).union(1))
+
+
+def test_issue_16469():
+    x = Symbol("x", real=True)
+    f = abs(x)
+    assert function_range(f, x, S.Reals) == Interval(0, oo, False, True)

sympy/calculus/util.py

@@ -10,6 +10,7 @@ from sympy.logic.boolalg import And
 from sympy.polys.rationaltools import together
 from sympy.sets.sets import (Interval, Intersection, FiniteSet, Union,
                              Complement, EmptySet)
+from sympy.sets.fancysets import ImageSet
 from sympy.simplify.radsimp import denom
 from sympy.solvers.inequalities import solve_univariate_inequality
 from sympy.utilities import filldedent
@@ -161,7 +162,7 @@ def function_range(f, symbol, domain):
         return S.EmptySet
 
     period = periodicity(f, symbol)
-    if period is S.Zero:
+    if period == S.Zero:
         # the expression is constant wrt symbol
         return FiniteSet(f.expand())
 
@@ -211,7 +212,11 @@ def function_range(f, symbol, domain):
             solution = solveset(f.diff(symbol), symbol, interval)
 
             if not iterable(solution):
-                raise NotImplementedError('Unable to find critical points for {}'.format(f))
+                raise NotImplementedError(
+                        'Unable to find critical points for {}'.format(f))
+            if isinstance(solution, ImageSet):
+                raise NotImplementedError(
+                        'Infinite number of critical points for {}'.format(f))
 
             critical_points += solution
 
@@ -269,7 +274,7 @@ def not_empty_in(finset_intersection, *syms):
     >>> not_empty_in(FiniteSet(x/2).intersect(Interval(0, 1)), x)
     Interval(0, 2)
     >>> not_empty_in(FiniteSet(x, x**2).intersect(Interval(1, 2)), x)
-    Union(Interval(-sqrt(2), -1), Interval(1, 2))
+    Union(Interval(1, 2), Interval(-sqrt(2), -1))
     >>> not_empty_in(FiniteSet(x**2/(x + 2)).intersect(Interval(1, oo)), x)
     Union(Interval.Lopen(-2, -1), Interval(2, oo))
     """
@@ -280,8 +285,8 @@ def not_empty_in(finset_intersection, *syms):
     if len(syms) == 0:
         raise ValueError("One or more symbols must be given in syms.")
 
-    if finset_intersection.is_EmptySet:
-        return EmptySet()
+    if finset_intersection is S.EmptySet:
+        return S.EmptySet
 
     if isinstance(finset_intersection, Union):
         elm_in_sets = finset_intersection.args[0]
@@ -406,7 +411,6 @@ def periodicity(f, symbol, check=False):
     pi
     >>> periodicity(exp(x), x)
     """
-    from sympy.core.function import diff
     from sympy.core.mod import Mod
     from sympy.core.relational import Relational
     from sympy.functions.elementary.exponential import exp
@@ -415,7 +419,7 @@ def periodicity(f, symbol, check=False):
         TrigonometricFunction, sin, cos, csc, sec)
     from sympy.simplify.simplify import simplify
     from sympy.solvers.decompogen import decompogen
-    from sympy.polys.polytools import degree, lcm_list
+    from sympy.polys.polytools import degree
 
     temp = Dummy('x', real=True)
     f = f.subs(symbol, temp)
@@ -582,7 +586,8 @@ def _periodicity(args, symbol):
     if len(periods) > 1:
         return lcim(periods)
 
-    return periods[0]
+    if periods:
+        return periods[0]
 
 
 def lcim(numbers):
@@ -995,17 +1000,15 @@ class AccumulationBounds(AtomicExpr):
     calculations, use ``mpmath.iv`` instead.
     """
 
-    is_real = True
+    is_extended_real = True
 
     def __new__(cls, min, max):
 
         min = _sympify(min)
         max = _sympify(max)
 
-        inftys = [S.Infinity, S.NegativeInfinity]
-        # Only allow real intervals (use symbols with 'is_real=True').
-        if not (min.is_real or min in inftys) \
-           or not (max.is_real or max in inftys):
+        # Only allow real intervals (use symbols with 'is_extended_real=True').
+        if not min.is_extended_real or not max.is_extended_real:
             raise ValueError("Only real AccumulationBounds are supported")
 
         # Make sure that the created AccumBounds object will be valid.
@@ -1022,6 +1025,10 @@ class AccumulationBounds(AtomicExpr):
     # setting the operation priority
     _op_priority = 11.0
 
+    def _eval_is_real(self):
+        if self.min.is_real and self.max.is_real:
+            return True
+
     @property
     def min(self):
         """
@@ -1102,7 +1109,7 @@ class AccumulationBounds(AtomicExpr):
             if other is S.Infinity and self.min is S.NegativeInfinity or \
                     other is S.NegativeInfinity and self.max is S.Infinity:
                 return AccumBounds(-oo, oo)
-            elif other.is_real:
+            elif other.is_extended_real:
                 return AccumBounds(Add(self.min, other), Add(self.max, other))
             return Add(self, other, evaluate=False)
         return NotImplemented
@@ -1122,7 +1129,7 @@ class AccumulationBounds(AtomicExpr):
             if other is S.NegativeInfinity and self.min is S.NegativeInfinity or \
                     other is S.Infinity and self.max is S.Infinity:
                 return AccumBounds(-oo, oo)
-            elif other.is_real:
+            elif other.is_extended_real:
                 return AccumBounds(
                     Add(self.min, -other),
                     Add(self.max, -other))
@@ -1155,7 +1162,7 @@ class AccumulationBounds(AtomicExpr):
                     return AccumBounds(-oo, 0)
                 if self.max.is_zero:
                     return AccumBounds(0, oo)
-            if other.is_real:
+            if other.is_extended_real:
                 if other.is_zero:
                     if self == AccumBounds(-oo, oo):
                         return AccumBounds(-oo, oo)
@@ -1164,11 +1171,11 @@ class AccumulationBounds(AtomicExpr):
                     if self.min is S.NegativeInfinity:
                         return AccumBounds(-oo, 0)
                     return S.Zero
-                if other.is_positive:
+                if other.is_extended_positive:
                     return AccumBounds(
                         Mul(self.min, other),
                         Mul(self.max, other))
-                elif other.is_negative:
+                elif other.is_extended_negative:
                     return AccumBounds(
                         Mul(self.max, other),
                         Mul(self.min, other))
@@ -1193,33 +1200,33 @@ class AccumulationBounds(AtomicExpr):
                         return AccumBounds(-oo, 0)
                     return AccumBounds(-oo, oo)
 
-                if self.max.is_negative:
-                    if other.min.is_negative:
+                if self.max.is_extended_negative:
+                    if other.min.is_extended_negative:
                         if other.max.is_zero:
                             return AccumBounds(self.max / other.min, oo)
-                        if other.max.is_positive:
+                        if other.max.is_extended_positive:
                             # the actual answer is a Union of AccumBounds,
                             # Union(AccumBounds(-oo, self.max/other.max),
                             #       AccumBounds(self.max/other.min, oo))
                             return AccumBounds(-oo, oo)
 
-                    if other.min.is_zero and other.max.is_positive:
+                    if other.min.is_zero and other.max.is_extended_positive:
                         return AccumBounds(-oo, self.max / other.max)
 
-                if self.min.is_positive:
-                    if other.min.is_negative:
+                if self.min.is_extended_positive:
+                    if other.min.is_extended_negative:
                         if other.max.is_zero:
                             return AccumBounds(-oo, self.min / other.min)
-                        if other.max.is_positive:
+                        if other.max.is_extended_positive:
                             # the actual answer is a Union of AccumBounds,
                             # Union(AccumBounds(-oo, self.min/other.min),
                             #       AccumBounds(self.min/other.max, oo))
                             return AccumBounds(-oo, oo)
 
-                    if other.min.is_zero and other.max.is_positive:
+                    if other.min.is_zero and other.max.is_extended_positive:
                         return AccumBounds(self.min / other.max, oo)
 
-            elif other.is_real:
+            elif other.is_extended_real:
                 if other is S.Infinity or other is S.NegativeInfinity:
                     if self == AccumBounds(-oo, oo):
                         return AccumBounds(-oo, oo)
@@ -1227,9 +1234,9 @@ class AccumulationBounds(AtomicExpr):
                         return AccumBounds(Min(0, other), Max(0, other))
                     if self.min is S.NegativeInfinity:
                         return AccumBounds(Min(0, -other), Max(0, -other))
-                if other.is_positive:
+                if other.is_extended_positive:
                     return AccumBounds(self.min / other, self.max / other)
-                elif other.is_negative:
+                elif other.is_extended_negative:
                     return AccumBounds(self.max / other, self.min / other)
             return Mul(self, 1 / other, evaluate=False)
 
@@ -1240,19 +1247,19 @@ class AccumulationBounds(AtomicExpr):
     @_sympifyit('other', NotImplemented)
     def __rdiv__(self, other):
         if isinstance(other, Expr):
-            if other.is_real:
+            if other.is_extended_real:
                 if other.is_zero:
                     return S.Zero
                 if S.Zero in self:
-                    if self.min == S.Zero:
-                        if other.is_positive:
+                    if self.min.is_zero:
+                        if other.is_extended_positive:
                             return AccumBounds(Mul(other, 1 / self.max), oo)
-                        if other.is_negative:
+                        if other.is_extended_negative:
                             return AccumBounds(-oo, Mul(other, 1 / self.max))
-                    if self.max == S.Zero:
-                        if other.is_positive:
+                    if self.max.is_zero:
+                        if other.is_extended_positive:
                             return AccumBounds(-oo, Mul(other, 1 / self.min))
-                        if other.is_negative:
+                        if other.is_extended_negative:
                             return AccumBounds(Mul(other, 1 / self.min), oo)
                     return AccumBounds(-oo, oo)
                 else:
@@ -1269,13 +1276,13 @@ class AccumulationBounds(AtomicExpr):
         from sympy.functions.elementary.miscellaneous import real_root
         if isinstance(other, Expr):
             if other is S.Infinity:
-                if self.min.is_nonnegative:
+                if self.min.is_extended_nonnegative:
                     if self.max < 1:
                         return S.Zero
                     if self.min > 1:
                         return S.Infinity
                     return AccumBounds(0, oo)
-                elif self.max.is_negative:
+                elif self.max.is_extended_negative:
                     if self.min > -1:
                         return S.Zero
                     if self.max < -1:
@@ -1291,22 +1298,22 @@ class AccumulationBounds(AtomicExpr):
             if other is S.NegativeInfinity:
                 return (1 / self)**oo
 
-            if other.is_real and other.is_number:
+            if other.is_extended_real and other.is_number:
                 if other.is_zero:
                     return S.One
 
                 if other.is_Integer:
-                    if self.min.is_positive:
+                    if self.min.is_extended_positive:
                         return AccumBounds(
                             Min(self.min ** other, self.max ** other),
                             Max(self.min ** other, self.max ** other))
-                    elif self.max.is_negative:
+                    elif self.max.is_extended_negative:
                         return AccumBounds(
                             Min(self.max ** other, self.min ** other),
                             Max(self.max ** other, self.min ** other))
 
                     if other % 2 == 0:
-                        if other.is_negative:
+                        if other.is_extended_negative:
                             if self.min.is_zero:
                                 return AccumBounds(self.max**other, oo)
                             if self.max.is_zero:
@@ -1315,7 +1322,7 @@ class AccumulationBounds(AtomicExpr):
                         return AccumBounds(
                             S.Zero, Max(self.min**other, self.max**other))
                     else:
-                        if other.is_negative:
+                        if other.is_extended_negative:
                             if self.min.is_zero:
                                 return AccumBounds(self.max**other, oo)
                             if self.max.is_zero:
@@ -1324,10 +1331,10 @@ class AccumulationBounds(AtomicExpr):
                         return AccumBounds(self.min**other, self.max**other)
 
                 num, den = other.as_numer_denom()
-                if num == S(1):
+                if num == S.One:
                     if den % 2 == 0:
                         if S.Zero in self:
-                            if self.min.is_negative:
+                            if self.min.is_extended_negative:
                                 return AccumBounds(0, real_root(self.max, den))
                     return AccumBounds(real_root(self.min, den),
                                        real_root(self.max, den))
@@ -1338,9 +1345,9 @@ class AccumulationBounds(AtomicExpr):
         return NotImplemented
 
     def __abs__(self):
-        if self.max.is_negative:
+        if self.max.is_extended_negative:
             return self.__neg__()
-        elif self.min.is_negative:
+        elif self.min.is_extended_negative:
             return AccumBounds(S.Zero, Max(abs(self.min), self.max))
         else:
             return self
@@ -1371,8 +1378,7 @@ class AccumulationBounds(AtomicExpr):
                 return True
             if self.min >= other.max:
                 return False
-        elif not(other.is_real or other is S.Infinity or
-               other is S.NegativeInfinity):
+        elif not other.is_extended_real:
             raise TypeError(
                 "Invalid comparison of %s %s" %
                 (type(other), other))
@@ -1409,8 +1415,7 @@ class AccumulationBounds(AtomicExpr):
                 return True
             if self.min > other.max:
                 return False
-        elif not(other.is_real or other is S.Infinity or
-               other is S.NegativeInfinity):
+        elif not other.is_extended_real:
             raise TypeError(
                 "Invalid comparison of %s %s" %
                 (type(other), other))
@@ -1447,8 +1452,7 @@ class AccumulationBounds(AtomicExpr):
                 return True
             if self.max <= other.min:
                 return False
-        elif not(other.is_real or other is S.Infinity or
-               other is S.NegativeInfinity):
+        elif not other.is_extended_real:
             raise TypeError(
                 "Invalid comparison of %s %s" %
                 (type(other), other))
@@ -1485,8 +1489,7 @@ class AccumulationBounds(AtomicExpr):
                 return True
             if self.max < other.min:
                 return False
-        elif not(other.is_real or other is S.Infinity or
-               other is S.NegativeInfinity):
+        elif not other.is_extended_real:
             raise TypeError(
                 "Invalid comparison of %s %s" %
                 (type(other), other))