code: def bezier2polynomial(p, numpy_ordering=True, return_poly1d=False):
      
      if len(p) == 4:
          coeffs = (-p[0] + 3*(p[1] - p[2]) + p[3],
                    3*(p[0] - 2*p[1] + p[2]),
                    3*(p[1]-p[0]),
                    p[0])
      elif len(p) == 3:
          coeffs = (p[0] - 2*p[1] + p[2],
                    2*(p[1] - p[0]),
                    p[0])
      elif len(p) == 2:
          coeffs = (p[1]-p[0],
                    p[0])
      elif len(p) == 1:
          coeffs = p
      else:
          # https://en.wikipedia.org/wiki/Bezier_curve#Polynomial_form
          n = len(p) - 1
          coeffs = [fac(n)//fac(n-j) * sum(
              (-1)**(i+j) * p[i] / (fac(i) * fac(j-i)) for i in range(j+1))
              for j in range(n+1)]
          coeffs.reverse()
      if not numpy_ordering:
          coeffs = coeffs[::-1]  # can't use .reverse() as might be tuple
      if return_poly1d:
          return poly1d(coeffs)
      return coeffs
docs: Converts a tuple of Bezier control points to a tuple of coefficients
      of the expanded polynomial.
      return_poly1d : returns a numpy.poly1d object.  This makes computations
      of derivatives/anti-derivatives and many other operations quite quick.
      numpy_ordering : By default (to accommodate numpy) the coefficients will
      be output in reverse standard order.