More efficient way to calculate magnetic field using Biot-Savart

I am writing a program in python that is supposed to calculate the magnetic field along a conducting coil that is made up of a bunch of points, and the magnetic field is generated by other conducting coils which were made by rotating the points that make up the original conducting coil. I have tried doing this using Biot-Savart law (code below), but each coil is made up of ~300 points, and there are 10 coils in total, so my method is called ~810,000 times, but I think creating and subtracting 2 vectors and then crossing them with each iteration makes it take a very long time to run, so I was wondering if someone knows of a more efficient way to do Biot-Savart, or if I should try a different method. Any insight would be very much appreciated, thank you!

My Biot-Savart method:

def biotsav(l,lprev,pointpos,cur,mu0=1e-7):

#l is the point along the wire generating the magnetic field
#lprev is the previous point along the wire that generated a magnetic field
#pointpos is the position where the magnetic field will be measured at
#cur is the current running through the wire

#if the point causing the magnetic field is the point where the magnetic field is being measured at, return 0 otherwise it would return infinity
if l[0] == pointpos[0] and l[1] == pointpos[1] and l[2] == pointpos[2]:
return np.array([0,0,0])

dl = np.subtract(l,lprev) #vector going from previous point along wire to current point
r = np.subtract(pointpos,l) #vector going from point causing magnetic field to point where magnetic field is being measured
rmag = np.sqrt(r[0]**2 + r[1]**2 + r[2]**2) #magnitude of r
return mu0*cur*np.cross(dl,r)/(rmag**3) #biot savart equation


The loop where the method is called:

#array that will hold magnetic field values along each point on the coil
Bfield = []

#Caclculates the magnetic field at a point on the TF caused by all points on all coils
for i in range(numpts): #this for loop runs through all points along original coil

Btemp = [0,0,0] #placeholder array for updating magnetic field measured at a point along original coil

for j in range(len(allxs)): #this for loop goes through all coils created through symmetry

for k in range(len(allxs[j])): #this for loop goes through all points along coils created by symmetry
currentpt = [allxs[j][k],allys[j][k],allzs[j][k]]
prevpt = [allxs[j][k-1],allys[j][k-1],allzs[j][k-1]]
coilpt = [xlist[i],ylist[i],zlist[i]]