simplify besselk_levin (use con. prop)

Author: heltonmc

src/BesselFunctions/besselk.jl

@@ -578,15 +578,14 @@ end
 @generated function besselkx_levin(v, x::T, ::Val{N}) where {T <: FloatTypes, N}
     :(
         begin
-            l = let (out, term, fv2) = (one(typeof(x)), one(typeof(x)), 4*v^2)
-                @ntuple $N i -> begin
-                    offset = muladd(2, i, -1)
-                    term *= muladd(offset, -offset, fv2) / (8 * x * i)
-                    out += term
-                    invterm = inv(term)
-                    Vec{2, T}((out * invterm, invterm))
+            s = zero(T)
+            t = one(T)
+            l = @ntuple $N i -> begin
+                    s += t
+                    t *= (4*v^2 - (2i - 1)^2) / (8 * x * i)
+                    invterm = inv(t)
+                    Vec{2, T}((s * invterm, invterm))
                 end
-            end
             return levin_transform(l) * sqrt(π / 2x)
         end
     )
@@ -595,16 +594,20 @@ end
 @generated function besselkx_levin(v, x::Complex{T}, ::Val{N}) where {T <: FloatTypes, N}
     :(
         begin
-            l = let (out, term, fv2, invx) = (one(typeof(x)), one(typeof(x)), 4*v^2, inv(8*x))
-                @ntuple $N i -> begin
-                    offset = muladd(2, i, -1)
-                    term *= @fastmath muladd(offset, -offset, fv2) * invx / i
-                    out += term
-                    invterm = @fastmath inv(term)
-                    Vec{4, T}((reim(out * invterm)..., reim(invterm)...))
+            s = zero(T)
+            t = one(typeof(x))
+            t2 = t
+            a = @fastmath inv(8*x)
+            a2 = 8*x
+
+            l =  @ntuple $N i -> begin
+                    s += t
+                    b = (4*v^2 - (2i - 1)^2) / i
+                    t *= a * b
+                    t2 *= a2 / b
+                    Vec{4, T}((reim(s * t2)..., reim(t2)...))
                 end
-            end
             return levin_transform(l) * sqrt(π / 2x)
         end
     )
-end
+end