code: def lagrange_interpolate(x, y, precision=250, **kwargs):
      
      n = len(x) - 1
      delta_x = [x2 - x1 for x1, x2 in zip(x, x[1:])]
      for i in range(n + 1):
          yield x[i], y[i]
          if i == n or delta_x[i] == 0:
              continue
  
          for s in range(1, precision):
              X = x[i] + s * delta_x[i] / precision
              s = 0
              for k in range(n + 1):
                  p = 1
                  for m in range(n + 1):
                      if m == k:
                          continue
                      if x[k] - x[m]:
                          p *= (X - x[m]) / (x[k] - x[m])
                  s += y[k] * p
              yield X, s
docs: Interpolate x, y using Lagrange polynomials
      https://en.wikipedia.org/wiki/Lagrange_polynomial