Merge pull request #60 from JuliaMath/mutate

Add mutating sequence

Author: heltonmc

src/besseli.jl

@@ -169,6 +169,8 @@ besseli(nu, x::Real) = _besseli(nu, float(x))
 
 _besseli(nu::Union{Int16, Float16}, x::Union{Int16, Float16}) = Float16(_besseli(Float32(nu), Float32(x)))
 
+_besseli(nu::AbstractRange, x::T) where T = besseli!(zeros(T, length(nu)), nu, x)
+
 function _besseli(nu::T, x::T) where T <: Union{Float32, Float64}
     isinteger(nu) && return _besseli(Int(nu), x)
     ~isfinite(x) && return x
@@ -206,11 +208,13 @@ function _besseli(nu::Integer, x::T) where T <: Union{Float32, Float64}
     end
 end
 
-function _besseli(nu::AbstractRange, x::T) where T
+besseli!(out::DenseVector, nu::AbstractRange, x) = _besseli!(out, nu, float(x))
+
+function _besseli!(out::DenseVector{T}, nu::AbstractRange, x::T) where T
     (nu[1] >= 0 && step(nu) == 1) || throw(ArgumentError("nu must be >= 0 with step(nu)=1"))
-    len = length(nu)
-    isone(len) && return [besseli(nu[1], x)]
-    out = zeros(T, len)
+    len = length(out)
+    !isequal(len, length(nu)) && throw(ArgumentError("out and nu must have the same length"))
+
     k = len
     inu = zero(T)
     while abs(inu) < floatmin(T)
@@ -225,8 +229,8 @@ function _besseli(nu::AbstractRange, x::T) where T
     end
     if k > 1
         out[k] = _besseli(nu[k], x)
-        tmp = @view out[1:k+1]
-        out[1:k+1] = besselk_down_recurrence!(tmp, x, nu[1:k+1])
+        tmp = @view out[begin:k+1]
+        besselk_down_recurrence!(tmp, x, nu[begin:k+1])
         return out
     else
         return out

src/besselj.jl

@@ -224,6 +224,8 @@ besselj(nu, x::Real) = _besselj(nu, float(x))
 
 _besselj(nu::Union{Int16, Float16}, x::Union{Int16, Float16}) = Float16(_besselj(Float32(nu), Float32(x)))
 
+_besselj(nu::AbstractRange, x::T) where T = besselj!(zeros(T, length(nu)), nu, x)
+
 function _besselj(nu::T, x::T) where T <: Union{Float32, Float64}
     isinteger(nu) && return _besselj(Int(nu), x)
     abs_nu = abs(nu)
@@ -269,14 +271,17 @@ function _besselj(nu::Integer, x::T) where T <: Union{Float32, Float64}
     end
 end
 
-function _besselj(nu::AbstractRange, x::T) where T
-    (nu[1] >= 0 && step(nu) == 1) || throw(ArgumentError("nu must be >= 0 with step(nu)=1"))
-    len = length(nu)
-    isone(len) && return [besselj(nu[1], x)]
+besselj!(out::DenseVector, nu::AbstractRange, x) = _besselj!(out, nu, float(x))
+
+function _besselj!(out::DenseVector{T}, nu::AbstractVector, x::T) where T <: Union{Float32, Float64}
+    (nu[begin] >= 0 && step(nu) == 1) || throw(ArgumentError("nu must be >= 0 with step(nu)=1"))
+    len = length(out)
+    !isequal(len, length(nu)) && throw(ArgumentError("out and nu must have the same length"))
+
+    isone(len) && return [besselj(nu[begin], x)]
 
-    out = zeros(T, len)
     if nu[end] < x
-        out[1], out[2] = _besselj(nu[1], x), _besselj(nu[2], x)
+        out[begin], out[begin + 1] = _besselj(nu[begin], x), _besselj(nu[begin + 1], x)
         return besselj_up_recurrence!(out, x, nu)
     else
         k = len
@@ -293,8 +298,8 @@ function _besselj(nu::AbstractRange, x::T) where T
         end
         if k > 1
             out[k] = _besselj(nu[k], x)
-            tmp = @view out[1:k+1]
-            besselj_down_recurrence!(tmp, x, nu[1:k+1])
+            tmp = @view out[begin:k+1]
+            besselj_down_recurrence!(tmp, x, nu[begin:k+1])
             return out
         else
             return out

src/besselk.jl

@@ -191,6 +191,8 @@ besselk(nu, x::Real) = _besselk(nu, float(x))
 
 _besselk(nu::Union{Int16, Float16}, x::Union{Int16, Float16}) = Float16(_besselk(Float32(nu), Float32(x)))
 
+_besselk(nu::AbstractRange, x::T) where T = besselk!(zeros(T, length(nu)), nu, x)
+
 function _besselk(nu::T, x::T) where T <: Union{Float32, Float64}
     isinteger(nu) && return _besselk(Int(nu), x)
     abs_nu = abs(nu)
@@ -216,11 +218,14 @@ function _besselk(nu::Integer, x::T) where T <: Union{Float32, Float64}
     end
 end
 
-function _besselk(nu::AbstractRange, x::T) where T
+besselk!(out::DenseVector, nu::AbstractRange, x) = _besselk!(out, nu, float(x))
+
+function _besselk!(out::DenseVector{T}, nu::AbstractRange, x::T) where T
     (nu[1] >= 0 && step(nu) == 1) || throw(ArgumentError("nu must be >= 0 with step(nu)=1"))
-    len = length(nu)
+    len = length(out)
+    !isequal(len, length(nu)) && throw(ArgumentError("out and nu must have the same length"))
     isone(len) && return [besselk(nu[1], x)]
-    out = zeros(T, len)
+
     k = 1
     knu = zero(T)
     while abs(knu) < floatmin(T)
@@ -236,7 +241,7 @@ function _besselk(nu::AbstractRange, x::T) where T
     if k < len
         out[k] = _besselk(nu[k], x)
         tmp = @view out[k-1:end]
-        out[k-1:end] = besselk_up_recurrence!(tmp, x, nu[k-1:end])
+        besselk_up_recurrence!(tmp, x, nu[k-1:end])
         return out
     else
         return out

src/bessely.jl

@@ -261,6 +261,8 @@ bessely(nu, x::Real) = _bessely(nu, float(x))
 
 _bessely(nu::Union{Int16, Float16}, x::Union{Int16, Float16}) = Float16(_bessely(Float32(nu), Float32(x)))
 
+_bessely(nu::AbstractRange, x::T) where T = bessely!(zeros(T, length(nu)), nu, x)
+
 function _bessely(nu::T, x::T) where T <: Union{Float32, Float64}
     isnan(nu) || isnan(x) && return NaN
     isinteger(nu) && return _bessely(Int(nu), x)
@@ -310,10 +312,12 @@ function _bessely(nu::Integer, x::T) where T <: Union{Float32, Float64}
     end
 end
 
-function _bessely(nu::AbstractRange, x::T) where T
-    (nu[1] >= 0 && step(nu) == 1) || throw(ArgumentError("nu must be >= 0 with step(nu)=1"))
-    out = Vector{T}(undef, length(nu))
-    out[1], out[2] = _bessely(nu[1], x), _bessely(nu[2], x)
+bessely!(out::DenseVector, nu::AbstractRange, x) = _bessely!(out, nu, float(x))
+
+function _bessely!(out::DenseVector{T}, nu::AbstractRange, x::T) where T
+    (nu[begin] >= 0 && step(nu) == 1) || throw(ArgumentError("nu must be >= 0 with step(nu)=1"))
+    !isequal(length(out), length(nu)) && throw(ArgumentError("out and nu must have the same length"))
+    out[begin], out[begin + 1] = _bessely(nu[begin], x), _bessely(nu[begin + 1], x)
     return besselj_up_recurrence!(out, x, nu)
 end
 

test/besseli_test.jl

@@ -109,6 +109,7 @@ end
 # test nu_range
 @test besseli(0:250, 2.0) ≈ SpecialFunctions.besseli.(0:250, 2.0) rtol=1e-13
 @test besseli(0.5:1:10.5, 2.0) ≈ SpecialFunctions.besseli.(0.5:1:10.5, 2.0) rtol=1e-13
+@test Bessels.besseli!(zeros(Float64, 10), 1:10, 1.0) ≈ besseli(1:10, 1.0)
 
 ### need to fix method ambiguities for other functions ###### 
 

test/besselj_test.jl

@@ -124,6 +124,7 @@ end
 @test besselj(0:95, 100.0) ≈ SpecialFunctions.besselj.(0:95, 100.0) rtol=1e-11
 @test besselj(0.5:1:150.5, 2.0) ≈ SpecialFunctions.besselj.(0.5:1:150.5, 2.0) rtol=1e-11
 @test besselj(0.5:1:10.5, 40.0) ≈ SpecialFunctions.besselj.(0.5:1:10.5, 40.0) rtol=1e-11
+@test Bessels.besselj!(zeros(Float64, 10), 1:10, 1.0) ≈ besselj(1:10, 1.0)
 
 # test Float16 and Float32
 @test besselj(Int16(10), Float16(1.0)) isa Float16

test/besselk_test.jl

@@ -130,6 +130,7 @@ end
 @test besselk(0:50, 2.0) ≈ SpecialFunctions.besselk.(0:50, 2.0) rtol=1e-13
 @test besselk(0.5:1:10.5, 12.0) ≈ SpecialFunctions.besselk.(0.5:1:10.5, 12.0) rtol=1e-13
 @test besselk(1:700, 800.0) ≈ SpecialFunctions.besselk.(1:700, 800.0)
+@test Bessels.besselk!(zeros(Float64, 10), 1:10, 1.0) ≈ besselk(1:10, 1.0)
 
 # test Float16
 @test besselk(Int16(10), Float16(1.0)) isa Float16

test/bessely_test.jl

@@ -88,7 +88,7 @@ end
 
 # test decimal orders
 # SpecialFunctions.jl can give errors over 1e-12 so need to soften tolerance to match
-# need to switch tests over to ArbNumerics.jl for better precision tests 
+# need to switch tests over to ArbNumerics.jl for better precision tests
 x = [0.05, 0.1, 0.2, 0.25, 0.3, 0.4, 0.5,0.55,  0.6,0.65,  0.7, 0.75, 0.8, 0.85, 0.9, 0.92, 0.95, 0.97, 0.99, 1.0, 1.01, 1.05, 1.08, 1.1, 1.2, 1.4, 1.5, 1.6, 1.8, 2.0, 2.5, 3.0, 4.0, 4.5, 4.99, 5.1]
 nu = [0.1, 0.4567, 0.8123, 1.5, 2.5, 4.1234, 6.8, 12.3, 18.9, 28.2345, 38.1235, 51.23, 72.23435, 80.5, 98.5, 104.2]
 for v in nu, xx in x
@@ -104,6 +104,7 @@ end
 @test bessely(0:50, 100.0) ≈ SpecialFunctions.bessely.(0:50, 100.0) rtol=1e-11
 @test bessely(0.5:1:10.5, 2.0) ≈ SpecialFunctions.bessely.(0.5:1:10.5, 2.0) rtol=1e-11
 @test bessely(0.5:1:10.5, 40.0) ≈ SpecialFunctions.bessely.(0.5:1:10.5, 40.0) rtol=1e-11
+@test Bessels.bessely!(zeros(Float64, 10), 1:10, 1.0) ≈ bessely(1:10, 1.0)
 
 # test Float16
 @test bessely(Int16(10), Float16(1.0)) isa Float16