code: def ginterp(coeff_matrix, points):
      
      # Dimension
      dim = len(points[0])
  
      # Number of data points
      num_points = len(points)
  
      # Solve system of linear equations
      matrix_l, matrix_u = linalg.lu_decomposition(coeff_matrix)
      ctrlpts = [[0.0 for _ in range(dim)] for _ in range(num_points)]
      for i in range(dim):
          b = [pt[i] for pt in points]
          y = linalg.forward_substitution(matrix_l, b)
          x = linalg.backward_substitution(matrix_u, y)
          for j in range(num_points):
              ctrlpts[j][i] = x[j]
  
      # Return control points
      return ctrlpts
docs: Applies global interpolation to the set of data points to find control points.
  
      :param coeff_matrix: coefficient matrix
      :type coeff_matrix: list, tuple
      :param points: data points
      :type points: list, tuple
      :return: control points
      :rtype: list