## Issue

So, I need help minimizing the time it takes to run the code with large numbers of data only by using NumPy. I think the for loops made my code inefficient.. But I do not know how to make the for loop into a list comprehension, which might help it run faster..

```
def lagrange(p,node,n,x):
m=[]
#base lagrange polynomial
for i in range(n):
for j in range(p+1):
L=1
for k in range(p+1):
if k!=j:
L= L*(x[i] - node[k])/(node[j] - node[k])
m.append(L)
lagrange= np.array(m).reshape(n,p+1)
return lagrange
def interpolant(a,b,p,n,x,f):
m=[]
node=np.linspace(a,b,p+1)
for j in range(n):
polynomial=0
for i in range(p+1):
polynomial += f(node[i]) * lagrange(p,node,n,x)
m.append(polynomial)
interpolant = np.array(inter)
return interpolant
```

## Solution

It appears the value of `lagrange_poly(...)`

is **recomputed n*(p+1) times for no reason** which is very very expensive! You can compute it once before the loop, store it in a variable and reuse the variable later.

Here is the fixed code:

```
def uniform_poly_interpolation(a,b,p,n,x,f,produce_fig):
inter=[]
xhat=np.linspace(a,b,p+1)
#use for loop to iterate interpolant.
mat = lagrange_poly(p,xhat,n,x,1e-10)[0]
for j in range(n):
po=0
for i in range(p+1):
po += f(xhat[i]) * mat[i,j]
inter.append(po)
interpolant = np.array(inter)
return interpolant
```

This should be much much faster.

Moreover, the execution is slow because **accessing scalar values of Numpy arrays from CPython is very slow**. Numpy is designed to work with array and not to extract scalar values in loops. Additionally, the loop CPython interpreter are relatively slow. You can solve this problem efficiently with **Numba** that compile your code to a very fast native code using a JIT-compiler.

Here is the Numba code:

```
import numba as nb
@nb.njit
def lagrange_poly(p, xhat, n, x, tol):
error_flag = 0
er = 1
lagrange_matrix = np.empty((n, p+1), dtype=np.float64)
for l in range(p):
if abs(xhat[l] - xhat[l+1]) < tol:
error_flag = er
# Base lagrange polynomial
for i in range(n):
for j in range(p+1):
L = 1.0
for k in range(p+1):
if k!=j:
L = L * (x[i] - xhat[k]) / (xhat[j] - xhat[k])
lagrange_matrix[i, j] = L
return lagrange_matrix, error_flag
```

Overall, this should be *several order of magnitude faster*.

Answered By – Jérôme Richard

Answer Checked By – Katrina (BugsFixing Volunteer)