mpmath · skirpichev · Apr 10, 2024 · Apr 6, 2024
diff --git a/mpmath/calculus/approximation.py b/mpmath/calculus/approximation.py
@@ -1,5 +1,8 @@
+import warnings
+
 from .calculus import defun
 
+
 #----------------------------------------------------------------------------#
 #                              Approximation methods                         #
 #----------------------------------------------------------------------------#
@@ -35,7 +38,7 @@ def chebT(ctx, a=1, b=0):
         Ta, Tb = Tmp, Ta
 
 @defun
-def chebyfit(ctx, f, interval, N, error=False):
+def chebyfit(ctx, f, interval, N, error=False, asc=None):
     r"""
     Computes a polynomial of degree `N-1` that approximates the
     given function `f` on the interval `[a, b]`. With ``error=True``,
@@ -54,30 +57,33 @@ def chebyfit(ctx, f, interval, N, error=False):
     example, it could be used to turn a slow mpmath function
     into a fast machine-precision version of the same.)
 
+    If *asc=False*, descending order of coefficients is used (the term
+    of largest degree - first).
+
     **Examples**
 
     Here we use :func:`~mpmath.chebyfit` to generate a low-degree approximation
     of `f(x) = \cos(x)`, valid on the interval `[1, 2]`::
 
         >>> from mpmath import mp, chebyfit, cos, nprint, polyval
         >>> mp.pretty = True
-        >>> poly, err = chebyfit(cos, [1, 2], 5, error=True)
+        >>> poly, err = chebyfit(cos, [1, 2], 5, error=True, asc=True)
         >>> nprint(poly)
-        [0.00291682, 0.146166, -0.732491, 0.174141, 0.949553]
+        [0.949553, 0.174141, -0.732491, 0.146166, 0.00291682]
         >>> nprint(err, 12)
         1.61351758081e-5
 
     The polynomial can be evaluated using ``polyval``::
 
-        >>> nprint(polyval(poly, 1.6), 12)
+        >>> nprint(polyval(poly, 1.6, asc=True), 12)
         -0.0291858904138
         >>> nprint(cos(1.6), 12)
         -0.0291995223013
 
     Sampling the true error at 1000 points shows that the error
     estimate generated by ``chebyfit`` is remarkably good::
 
-        >>> error = lambda x: abs(cos(x) - polyval(poly, x))
+        >>> error = lambda x: abs(cos(x) - polyval(poly, x, asc=True))
         >>> nprint(max([error(1+n/1000.) for n in range(1000)]), 12)
         1.61349954245e-5
 
@@ -122,18 +128,23 @@ def chebyfit(ctx, f, interval, N, error=False):
         for (k, Tk) in zip(range(N), T):
             for i in range(len(Tk)):
                 d[i] += c[k]*Tk[i]
-        d = d[::-1]
         # Estimate maximum error
         err = ctx.zero
         for k in range(N):
             x = ctx.cos(ctx.pi*k/N) * (b-a)*h + (b+a)*h
-            err = max(err, abs(f(x) - ctx.polyval(d, x)))
+            err = max(err, abs(f(x) - ctx.polyval(d, x, asc=True)))
     finally:
         ctx.prec = orig
+    if asc is None:
+        warnings.warn("Descending (wrt powers) order of polynomial "
+                      "coefficients is deprecated, please adapt you "
+                      "code to use ascending order, asc=True.",
+                      DeprecationWarning)
+        asc = False
     if error:
-        return d, +err
+        return d if asc else d[::-1], +err
     else:
-        return d
+        return d if asc else d[::-1]
 
 @defun
 def fourier(ctx, f, interval, N):

diff --git a/mpmath/calculus/differentiation.py b/mpmath/calculus/differentiation.py
@@ -1,5 +1,6 @@
 from .calculus import defun
 
+
 try:
     iteritems = dict.iteritems
 except AttributeError:
@@ -559,12 +560,11 @@ def taylor(ctx, f, x, n, **options):
     and supported keyword options.
 
     Note that to evaluate the Taylor polynomial as an approximation
-    of `f`, e.g. with :func:`~mpmath.polyval`, the coefficients must be reversed,
-    and the point of the Taylor expansion must be subtracted from
+    of `f`, the point of the Taylor expansion must be subtracted from
     the argument:
 
         >>> p = taylor(exp, 2.0, 10)
-        >>> polyval(p[::-1], 2.5 - 2.0)
+        >>> polyval(p, 2.5 - 2.0, asc=True)
         12.1824939606092
         >>> exp(2.5)
         12.1824939607035
@@ -608,7 +608,7 @@ def pade(ctx, a, L, M):
         >>> a = taylor(f, 0, 6)
         >>> p, q = pade(a, 3, 3)
         >>> x = 10
-        >>> polyval(p[::-1], x)/polyval(q[::-1], x)
+        >>> polyval(p, x, asc=True)/polyval(q, x, asc=True)
         1.38169105566806
         >>> f(x)
         1.38169855941551

diff --git a/mpmath/calculus/odes.py b/mpmath/calculus/odes.py
@@ -1,5 +1,6 @@
 from bisect import bisect
 
+
 class ODEMethods:
     pass
 
@@ -246,7 +247,7 @@ def odefun(ctx, F, x0, y0, tol=None, degree=None, method='taylor', verbose=False
     series_data = [(ser, x0, xb)]
     # We will be working with vectors of Taylor series
     def mpolyval(ser, a):
-        return [ctx.polyval(s[::-1], a) for s in ser]
+        return [ctx.polyval(s, a, asc=True) for s in ser]
     # Find nearest expansion point; compute if necessary
     def get_series(x):
         if x < x0:

diff --git a/mpmath/calculus/polynomials.py b/mpmath/calculus/polynomials.py
@@ -1,39 +1,52 @@
+import warnings
+
 from .calculus import defun
 
+
 #----------------------------------------------------------------------------#
 #                                Polynomials                                 #
 #----------------------------------------------------------------------------#
 
 # XXX: extra precision
 @defun
-def polyval(ctx, coeffs, x, derivative=False):
+def polyval(ctx, coeffs, x, derivative=False, asc=None):
     r"""
-    Given coefficients `[c_n, \ldots, c_2, c_1, c_0]` and a number `x`,
+    Given coefficients `[c_0, c_1, c_2, \ldots, c_n]` and a number `x`,
     :func:`~mpmath.polyval` evaluates the polynomial
 
     .. math ::
 
-        P(x) = c_n x^n + \ldots + c_2 x^2 + c_1 x + c_0.
+        P(x) = c_0 + c_1 x + c_2 x^2 \ldots c_n x^n
 
     If *derivative=True* is set, :func:`~mpmath.polyval` simultaneously
     evaluates `P(x)` with the derivative, `P'(x)`, and returns the
     tuple `(P(x), P'(x))`.
 
         >>> from mpmath import mp, polyval
         >>> mp.pretty = True
-        >>> polyval([3, 0, 2], 0.5)
+        >>> polyval([2, 0, 3], 0.5, asc=True)
         2.75
-        >>> polyval([3, 0, 2], 0.5, derivative=True)
+        >>> polyval([2, 0, 3], 0.5, derivative=True, asc=True)
         (2.75, 3.0)
 
+    If *asc=False*, descending order of coefficients is used (the term
+    of largest degree - first).
+
     The coefficients and the evaluation point may be any combination
     of real or complex numbers.
     """
     if not coeffs:
         return ctx.zero
-    p = ctx.convert(coeffs[0])
+    if asc is None:
+        warnings.warn("Descending (wrt powers) order of polynomial "
+                      "coefficients is deprecated, please adapt you "
+                      "code to use ascending order, asc=True.",
+                      DeprecationWarning)
+        asc = False
+        coeffs = coeffs[::-1]
+    p = ctx.convert(coeffs[-1])
     q = ctx.zero
-    for c in coeffs[1:]:
+    for c in reversed(coeffs[:-1]):
         if derivative:
             q = p + x*q
         p = c + x*p
@@ -44,7 +57,7 @@ def polyval(ctx, coeffs, x, derivative=False):
 
 @defun
 def polyroots(ctx, coeffs, maxsteps=50, cleanup=True, extraprec=10,
-        error=False, roots_init=None):
+              error=False, roots_init=None, asc=None):
     """
     Computes all roots (real or complex) of a given polynomial.
 
@@ -56,19 +69,22 @@ def polyroots(ctx, coeffs, maxsteps=50, cleanup=True, extraprec=10,
     With *error=True*, :func:`~mpmath.polyroots` returns a tuple *(roots, err)*
     where *err* is an estimate of the maximum error among the computed roots.
 
+    If *asc=False*, descending order of coefficients is used (the term
+    of largest degree - first).
+
     **Examples**
 
     Finding the three real roots of `x^3 - x^2 - 14x + 24`::
 
         >>> from mpmath import mp, polyroots, nprint, sqrt, polyval
         >>> mp.pretty = True
-        >>> nprint(polyroots([1,-1,-14,24]), 4)
+        >>> nprint(polyroots([24,-14,-1,1],asc=True), 4)
         [-4.0, 2.0, 3.0]
 
     Finding the two complex conjugate roots of `4x^2 + 3x + 2`, with an
     error estimate::
 
-        >>> roots, err = polyroots([4,3,2], error=True)
+        >>> roots, err = polyroots([2,3,4], error=True, asc=True)
         >>> for r in roots:
         ...     print(r)
         ...
@@ -78,16 +94,16 @@ def polyroots(ctx, coeffs, maxsteps=50, cleanup=True, extraprec=10,
         >>> err
         2.22044604925031e-16
         >>>
-        >>> polyval([4,3,2], roots[0])
+        >>> polyval([2,3,4], roots[0], asc=True)
         (2.22044604925031e-16 + 0.0j)
-        >>> polyval([4,3,2], roots[1])
+        >>> polyval([2,3,4], roots[1], asc=True)
         (2.22044604925031e-16 + 0.0j)
 
     The following example computes all the 5th roots of unity; that is,
     the roots of `x^5 - 1`::
 
         >>> mp.dps = 20
-        >>> for r in polyroots([1, 0, 0, 0, 0, -1]):
+        >>> for r in polyroots([-1, 0, 0, 0, 0, 1], asc=True):
         ...     print(r)
         ...
         1.0
@@ -116,7 +132,7 @@ def polyroots(ctx, coeffs, maxsteps=50, cleanup=True, extraprec=10,
     typically compute all roots of an arbitrary polynomial to high precision::
 
         >>> mp.dps = 60
-        >>> for r in polyroots([1, 0, -10, 0, 1]):
+        >>> for r in polyroots([1, 0, -10, 0, 1], asc=True):
         ...     print(r)
         ...
         -3.14626436994197234232913506571557044551247712918732870123249
@@ -156,17 +172,25 @@ def polyroots(ctx, coeffs, maxsteps=50, cleanup=True, extraprec=10,
         # Constant polynomial with no roots
         return []
 
+    if asc is None:
+        warnings.warn("Descending (wrt powers) order of polynomial "
+                      "coefficients is deprecated, please adapt you "
+                      "code to use ascending order, asc=True.",
+                      DeprecationWarning)
+        asc = False
+        coeffs = coeffs[::-1]
+
     orig = ctx.prec
     tol = +ctx.eps
     with ctx.extraprec(extraprec):
         deg = len(coeffs) - 1
         # Must be monic
-        lead = ctx.convert(coeffs[0])
+        lead = ctx.convert(coeffs[-1])
         if lead == 1:
             coeffs = [ctx.convert(c) for c in coeffs]
         else:
             coeffs = [c/lead for c in coeffs]
-        f = lambda x: ctx.polyval(coeffs, x)
+        f = lambda x: ctx.polyval(coeffs, x, asc=True)
         if roots_init is None:
             roots = [ctx.mpc((0.4+0.9j)**n) for n in range(deg)]
         else:

diff --git a/mpmath/function_docs.py b/mpmath/function_docs.py
@@ -2141,7 +2141,7 @@
 
     >>> nsum(lambda k: 1/(k**2-2*k+3), [0, inf])
     1.694361433907061256154665
-    >>> nprint(polyroots([1,-2,3]))
+    >>> nprint(polyroots([3,-2,1], asc=True))
     [(1.0 - 1.41421j), (1.0 + 1.41421j)]
     >>> r1 = 1-sqrt(2)*j
     >>> r2 = r1.conjugate()
@@ -5550,7 +5550,7 @@
 on the interval `[-1, 1]`::
 
     >>> for n in range(5):
-    ...     nprint(polyroots(taylor(lambda x: legendre(n, x), 0, n)[::-1]))
+    ...     nprint(polyroots(taylor(lambda x: legendre(n, x), 0, n), asc=True))
     ...
     []
     [0.0]
@@ -7606,7 +7606,7 @@
 can be checked to agree with the list of primitive roots::
 
     >>> p = taylor(lambda x: cyclotomic(6,x), 0, 6)[:3]
-    >>> for r in polyroots(p[::-1]):
+    >>> for r in polyroots(p, asc=True):
     ...     print(r)
     ...
     (0.5 - 0.8660254037844386467637232j)