diff --git a/doc/release/1.12.0-notes.rst b/doc/release/1.12.0-notes.rst
index 07871c10ceef..6a6cfcfcb342 100644
--- a/doc/release/1.12.0-notes.rst
+++ b/doc/release/1.12.0-notes.rst
@@ -128,7 +128,7 @@ file that will remain empty (bar a docstring) in the standard numpy source,
 but that can be overwritten by people making binary distributions of numpy.
 
 New nanfunctions ``nancumsum`` and ``nancumprod`` added
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 Nanfunctions ``nancumsum`` and ``nancumprod`` have been added to
 compute ``cumsum`` and ``cumprod`` by ignoring nans.
 
@@ -137,6 +137,13 @@ compute ``cumsum`` and ``cumprod`` by ignoring nans.
 ``np.lib.interp(x, xp, fp)`` now allows the interpolated array ``fp``
 to be complex and will interpolate at ``complex128`` precision.
 
+New polynomial evaluation function ``polyvalfromroots`` added
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+The new function ``polyvalfromroots`` evaluates a polynomial at given points
+from the roots of the polynomial. This is useful for higher order polynomials,
+where expansion into polynomial coefficients is inaccurate at machine
+precision.
+
 Improvements
 ============
 
diff --git a/doc/source/reference/routines.polynomials.polynomial.rst b/doc/source/reference/routines.polynomials.polynomial.rst
index 431856622bca..8194ca867fc1 100644
--- a/doc/source/reference/routines.polynomials.polynomial.rst
+++ b/doc/source/reference/routines.polynomials.polynomial.rst
@@ -32,6 +32,7 @@ Basics
    polygrid3d
    polyroots
    polyfromroots
+   polyvalfromroots
 
 Fitting
 -------
diff --git a/numpy/polynomial/polynomial.py b/numpy/polynomial/polynomial.py
index 5d05f599191e..23f1f6a36329 100644
--- a/numpy/polynomial/polynomial.py
+++ b/numpy/polynomial/polynomial.py
@@ -36,6 +36,7 @@
 --------------
 - `polyfromroots` -- create a polynomial with specified roots.
 - `polyroots` -- find the roots of a polynomial.
+- `polyvalfromroots` -- evalute a polynomial at given points from roots.
 - `polyvander` -- Vandermonde-like matrix for powers.
 - `polyvander2d` -- Vandermonde-like matrix for 2D power series.
 - `polyvander3d` -- Vandermonde-like matrix for 3D power series.
@@ -58,8 +59,8 @@
 __all__ = [
     'polyzero', 'polyone', 'polyx', 'polydomain', 'polyline', 'polyadd',
     'polysub', 'polymulx', 'polymul', 'polydiv', 'polypow', 'polyval',
-    'polyder', 'polyint', 'polyfromroots', 'polyvander', 'polyfit',
-    'polytrim', 'polyroots', 'Polynomial', 'polyval2d', 'polyval3d',
+    'polyvalfromroots', 'polyder', 'polyint', 'polyfromroots', 'polyvander',
+    'polyfit', 'polytrim', 'polyroots', 'Polynomial', 'polyval2d', 'polyval3d',
     'polygrid2d', 'polygrid3d', 'polyvander2d', 'polyvander3d']
 
 import warnings
@@ -780,6 +781,93 @@ def polyval(x, c, tensor=True):
     return c0
 
 
+def polyvalfromroots(x, r, tensor=True):
+    """
+    Evaluate a polynomial specified by its roots at points x.
+
+    If `r` is of length `N`, this function returns the value
+
+    .. math:: p(x) = \prod_{n=1}^{N} (x - r_n)
+
+    The parameter `x` is converted to an array only if it is a tuple or a
+    list, otherwise it is treated as a scalar. In either case, either `x`
+    or its elements must support multiplication and addition both with
+    themselves and with the elements of `r`.
+
+    If `r` is a 1-D array, then `p(x)` will have the same shape as `x`.  If `r`
+    is multidimensional, then the shape of the result depends on the value of
+    `tensor`. If `tensor is ``True`` the shape will be r.shape[1:] + x.shape;
+    that is, each polynomial is evaluated at every value of `x`. If `tensor` is
+    ``False``, the shape will be r.shape[1:]; that is, each polynomial is
+    evaluated only for the corresponding broadcast value of `x`. Note that
+    scalars have shape (,).
+
+    .. versionadded:: 1.12
+
+    Parameters
+    ----------
+    x : array_like, compatible object
+        If `x` is a list or tuple, it is converted to an ndarray, otherwise
+        it is left unchanged and treated as a scalar. In either case, `x`
+        or its elements must support addition and multiplication with
+        with themselves and with the elements of `r`.
+    r : array_like
+        Array of roots. If `r` is multidimensional the first index is the
+        root index, while the remaining indices enumerate multiple
+        polynomials. For instance, in the two dimensional case the roots
+        of each polynomial may be thought of as stored in the columns of `r`.
+    tensor : boolean, optional
+        If True, the shape of the roots array is extended with ones on the
+        right, one for each dimension of `x`. Scalars have dimension 0 for this
+        action. The result is that every column of coefficients in `r` is
+        evaluated for every element of `x`. If False, `x` is broadcast over the
+        columns of `r` for the evaluation.  This keyword is useful when `r` is
+        multidimensional. The default value is True.
+
+    Returns
+    -------
+    values : ndarray, compatible object
+        The shape of the returned array is described above.
+
+    See Also
+    --------
+    polyroots, polyfromroots, polyval
+
+    Examples
+    --------
+    >>> from numpy.polynomial.polynomial import polyvalfromroots
+    >>> polyvalfromroots(1, [1,2,3])
+    0.0
+    >>> a = np.arange(4).reshape(2,2)
+    >>> a
+    array([[0, 1],
+           [2, 3]])
+    >>> polyvalfromroots(a, [-1, 0, 1])
+    array([[ -0.,   0.],
+           [  6.,  24.]])
+    >>> r = np.arange(-2, 2).reshape(2,2) # multidimensional coefficients
+    >>> r # each column of r defines one polynomial
+    array([[-2, -1],
+           [ 0,  1]])
+    >>> b = [-2, 1]
+    >>> polyvalfromroots(b, r, tensor=True)
+    array([[-0.,  3.],
+           [ 3., 0.]])
+    >>> polyvalfromroots(b, r, tensor=False)
+    array([-0.,  0.])
+    """
+    r = np.array(r, ndmin=1, copy=0) + 0.0
+    if r.size == 0:
+        return np.ones_like(x)
+    if np.isscalar(x) or isinstance(x, (tuple, list)):
+        x = np.array(x, ndmin=1, copy=0) + 0.0
+    if isinstance(x, np.ndarray) and tensor:
+        if x.size == 0:
+            return x.reshape(r.shape[1:]+x.shape)
+        r = r.reshape(r.shape + (1,)*x.ndim)
+    return np.prod(x - r, axis=0)
+
+
 def polyval2d(x, y, c):
     """
     Evaluate a 2-D polynomial at points (x, y).
diff --git a/numpy/polynomial/tests/test_polynomial.py b/numpy/polynomial/tests/test_polynomial.py
index 0e6a2e8a006e..10c3a9c2ea46 100644
--- a/numpy/polynomial/tests/test_polynomial.py
+++ b/numpy/polynomial/tests/test_polynomial.py
@@ -136,6 +136,65 @@ def test_polyval(self):
             assert_equal(poly.polyval(x, [1, 0]).shape, dims)
             assert_equal(poly.polyval(x, [1, 0, 0]).shape, dims)
 
+    def test_polyvalfromroots(self):
+        #check empty input
+        assert_equal(poly.polyvalfromroots([], [1]).size, 0)
+        assert_(poly.polyvalfromroots([], [1]).shape == (0,))
+
+        #check empty input + multidimensional roots
+        assert_equal(poly.polyvalfromroots([], [[1] * 5]).size, 0)
+        assert_(poly.polyvalfromroots([], [[1] * 5]).shape == (5, 0))
+
+        #check scalar input
+        assert_equal(poly.polyvalfromroots(1, 1), 0)
+        assert_(poly.polyvalfromroots(1, np.ones((3, 3))).shape == (3, 1))
+
+        #check normal input)
+        x = np.linspace(-1, 1)
+        y = [x**i for i in range(5)]
+        for i in range(1, 5):
+            tgt = y[i]
+            res = poly.polyvalfromroots(x, [0]*i)
+            assert_almost_equal(res, tgt)
+        tgt = x*(x - 1)*(x + 1)
+        res = poly.polyvalfromroots(x, [-1, 0, 1])
+        assert_almost_equal(res, tgt)
+
+        #check that shape is preserved
+        for i in range(3):
+            dims = [2]*i
+            x = np.zeros(dims)
+            assert_equal(poly.polyvalfromroots(x, [1]).shape, dims)
+            assert_equal(poly.polyvalfromroots(x, [1, 0]).shape, dims)
+            assert_equal(poly.polyvalfromroots(x, [1, 0, 0]).shape, dims)
+
+        #check compatibility with factorization
+        ptest = [15, 2, -16, -2, 1]
+        r = poly.polyroots(ptest)
+        x = np.linspace(-1, 1)
+        assert_almost_equal(poly.polyval(x, ptest),
+                            poly.polyvalfromroots(x, r))
+
+        #check multidimensional arrays of roots and values
+        #check tensor=False
+        rshape = (3, 5)
+        x = np.arange(-3, 2)
+        r = np.random.randint(-5, 5, size=rshape)
+        res = poly.polyvalfromroots(x, r, tensor=False)
+        tgt = np.empty(r.shape[1:])
+        for ii in range(tgt.size):
+            tgt[ii] = poly.polyvalfromroots(x[ii], r[:, ii])
+        assert_equal(res, tgt)
+
+        #check tensor=True
+        x = np.vstack([x, 2*x])
+        res = poly.polyvalfromroots(x, r, tensor=True)
+        tgt = np.empty(r.shape[1:] + x.shape)
+        for ii in range(r.shape[1]):
+            for jj in range(x.shape[0]):
+                tgt[ii, jj, :] = poly.polyvalfromroots(x[jj], r[:, ii])
+        assert_equal(res, tgt)
+
     def test_polyval2d(self):
         x1, x2, x3 = self.x
         y1, y2, y3 = self.y