code: def chebyfit(ctx, f, interval, N, error=False):
      r"""
      Computes a polynomial of degree `N-1` that approximates the
      given function `f` on the interval `[a, b]`. With ``error=True``,
      :func:`~mpmath.chebyfit` also returns an accurate estimate of the
      maximum absolute error; that is, the maximum value of
      `|f(x) - P(x)|` for `x \in [a, b]`.
  
      :func:`~mpmath.chebyfit` uses the Chebyshev approximation formula,
      which gives a nearly optimal solution: that is, the maximum
      error of the approximating polynomial is very close to
      the smallest possible for any polynomial of the same degree.
  
      Chebyshev approximation is very useful if one needs repeated
      evaluation of an expensive function, such as function defined
      implicitly by an integral or a differential equation. (For
      example, it could be used to turn a slow mpmath function
      into a fast machine-precision version of the same.)
  
      **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.dps = 15; mp.pretty = True
          >>> poly, err = chebyfit(cos, [1, 2], 5, error=True)
          >>> nprint(poly)
          [0.00291682, 0.146166, -0.732491, 0.174141, 0.949553]
          >>> nprint(err, 12)
          1.61351758081e-5
  
      The polynomial can be evaluated using ``polyval``::
  
          >>> nprint(polyval(poly, 1.6), 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))
          >>> nprint(max([error(1+n/1000.) for n in range(1000)]), 12)
          1.61349954245e-5
  
      **Choice of degree**
  
      The degree `N` can be set arbitrarily high, to obtain an
      arbitrarily good approximation. As a rule of thumb, an
      `N`-term Chebyshev approximation is good to `N/(b-a)` decimal
      places on a unit interval (although this depends on how
      well-behaved `f` is). The cost grows accordingly: ``chebyfit``
      evaluates the function `(N^2)/2` times to compute the
      coefficients and an additional `N` times to estimate the error.
  
      **Possible issues**
  
      One should be careful to use a sufficiently high working
      precision both when calling ``chebyfit`` and when evaluating
      the resulting polynomial, as the polynomial is sometimes
      ill-conditioned. It is for example difficult to reach
      15-digit accuracy when evaluating the polynomial using
      machine precision floats, no matter the theoretical
      accuracy of the polynomial. (The option to return the
      coefficients in Chebyshev form should be made available
      in the future.)
  
      It is important to note the Chebyshev approximation works
      poorly if `f` is not smooth. A function containing singularities,
      rapid oscillation, etc can be approximated more effectively by
      multiplying it by a weight function that cancels out the
      nonsmooth features, or by dividing the interval into several
      segments.
      """
      a, b = ctx._as_points(interval)
      orig = ctx.prec
      try:
          ctx.prec = orig + int(N**0.5) + 20
          c = [chebcoeff(ctx, f, a, b, k, N) for k in range(N)]
          d = [ctx.zero] * N
          d[0] = -c[0] / 2
          h = ctx.mpf(0.5)
          T = chebT(
              ctx, ctx.mpf(2) / (b - a), ctx.mpf(-1) * (b + a) / (b - a)
          )
          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)))
      finally:
          ctx.prec = orig
      if error:
          return d, +err
      else:
          return d
  
imports: from ..libmp.backend import xrange
  from .calculus import defun
label_desc: Use image node_id 3 for calling a global function with example usage: chebyfit(ctx, f, interval, N, error) and returns: d,  unknown, d