Propagate inbounds in dnPl (#6)

* remove OffsetArray wrapper in collectPl!

* propagate inbounds if possible

Author: jishnub

Project.toml

@@ -1,6 +1,6 @@
 name = "LegendrePolynomials"
 uuid = "3db4a2ba-fc88-11e8-3e01-49c72059a882"
-version = "0.3.2"
+version = "0.3.3"
 
 [deps]
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"

src/LegendrePolynomials.jl

@@ -16,9 +16,6 @@ end
 
 assertnonnegative(l) = (l >= 0 || throw(ArgumentError("l must be >= 0, received " * string(l))))
 
-function checksize(arr, lmax)
-	maximum(axes(arr,1)) >= lmax || throw(ArgumentError("array is not large enough to store all values"))
-end
 function checklength(arr, minlength)
 	length(arr) >= minlength || throw(ArgumentError(
 		"array is not large enough to store all values, require a minimum length of " * string(minlength)))
@@ -168,6 +165,39 @@ function _checkvalues(x, l, n)
 	n >= 0 || throw(ArgumentError("order of derivative n must be >= 0"))
 end
 
+Base.@propagate_inbounds function _unsafednPl!(cache, x, l, n)
+	# unsafe, assumes 1-based indexing
+	checklength(cache, l - n + 1)
+	if n == l # may short-circuit this
+		cache[1] = doublefactorial(eltype(cache), 2l-1)
+	else
+		collectPl!(cache, x, lmax = l - n)
+
+		for ni in 1:n
+			# We denote the terms as P_ni_li
+
+			# li == ni
+			P_nim1_nim1 = cache[1]
+			P_ni_ni = dPl_recursion(eltype(cache), ni, ni, nothing, P_nim1_nim1, nothing, x)
+			cache[1] = P_ni_ni
+
+			# li == ni + 1
+			P_nim1_ni = cache[2]
+			P_ni_nip1 = dPl_recursion(eltype(cache), ni + 1, ni, P_ni_ni, P_nim1_ni, nothing, x)
+			cache[2] = P_ni_nip1
+
+			for li in ni+2:min(l, l - n + ni)
+				P_ni_lim2 = cache[li - ni - 1]
+				P_ni_lim1 = cache[li - ni]
+				P_nim1_lim1 = cache[li - ni + 1]
+				P_ni_li = dPl_recursion(eltype(cache), li, ni, P_ni_lim1, P_nim1_lim1, P_ni_lim2, x)
+				cache[li - ni + 1] = P_ni_li
+			end
+		end
+	end
+	nothing
+end
+
 """
 	dnPl(x, l::Integer, n::Integer, [cache::AbstractVector])
 
@@ -188,7 +218,7 @@ julia> dnPl(0.5, 4, 0) == Pl(0.5, 4) # zero-th order derivative == Pl(x)
 true
 ```
 """
-function dnPl(x, l::Integer, n::Integer, 
+Base.@propagate_inbounds function dnPl(x, l::Integer, n::Integer, 
 	A = begin
 		_checkvalues(x, l, n)
 		# do not allocate A if the value is trivially zero
@@ -200,42 +230,14 @@ function dnPl(x, l::Integer, n::Integer,
 	)
 
 	_checkvalues(x, l, n)
-	checklength(A, l - n + 1)
-
-	cache = OffsetArrays.no_offset_view(A)
-
 	# check if the value is trivially zero in case A is provided in the function call
 	if l < n
-		return zero(eltype(cache))
+		return zero(eltype(A))
 	end
-	
-	if n == l # may short-circuit this
-		cache[1] = doublefactorial(eltype(cache), 2l-1)
-	else
-		collectPl!(cache, x, lmax = l - n)
-
-		for ni in 1:n
-			# We denote the terms as P_ni_li
-
-			# li == ni
-			P_nim1_nim1 = cache[1]
-			P_ni_ni = dPl_recursion(eltype(cache), ni, ni, nothing, P_nim1_nim1, nothing, x)
-			cache[1] = P_ni_ni
 
-			# li == ni + 1
-			P_nim1_ni = cache[2]
-			P_ni_nip1 = dPl_recursion(eltype(cache), ni + 1, ni, P_ni_ni, P_nim1_ni, nothing, x)
-			cache[2] = P_ni_nip1
-
-			for li in ni+2:min(l, l - n + ni)
-				P_ni_lim2 = cache[li - ni - 1]
-				P_ni_lim1 = cache[li - ni]
-				P_nim1_lim1 = cache[li - ni + 1]
-				P_ni_li = dPl_recursion(eltype(cache), li, ni, P_ni_lim1, P_nim1_lim1, P_ni_lim2, x)
-				cache[li - ni + 1] = P_ni_li
-			end
-		end
-	end
+	cache = OffsetArrays.no_offset_view(A)
+	# function barrier, as no_offset_view may be type-unstable
+	_unsafednPl!(cache, x, l, n)
 
 	return cache[l - n + 1]
 end
@@ -271,13 +273,12 @@ julia> collectPl!(v, 0.5, lmax = 3) # only l from 0 to 3 are populated
 ```
 """
 function collectPl!(v::AbstractVector, x; lmax::Integer  = length(v) - 1)
-	v_0based = OffsetArray(v, OffsetArrays.Origin(0))
-	checksize(v_0based, lmax)
+	checklength(v, lmax + 1)
 
 	iter = LegendrePolynomialIterator(x, lmax)
 
 	for (l, Pl) in pairs(iter)
-		v_0based[l] = Pl
+		v[l + firstindex(v)] = Pl
 	end
 
 	v
@@ -363,7 +364,7 @@ function collectdnPl!(v, x; lmax::Integer, n::Integer)
 	# trivially zero for l < n
 	fill!((@view v[(0:n-1) .+ firstindex(v)]), zero(eltype(v)))
 	# populate the other elements
-	dnPl(x, lmax, n, @view v[(n:lmax) .+ firstindex(v)])
+	@inbounds dnPl(x, lmax, n, @view v[(n:lmax) .+ firstindex(v)])
 
 	v
 end