Implement select(k) using quickselect

Separate stats stuff into statistics.j
Implement median using select

j/linalg.j

@@ -237,18 +237,6 @@ function trace{T}(A::Matrix{T})
     return t
 end
 
-mean(V::Vector) = sum(V) / length(V)
-
-function std(V::Vector)
-    n = numel(V)
-    m = mean(V)
-    s = 0.0
-    for i=1:n
-        s += (V[i] - m)^2
-    end
-    return sqrt(s/(n-1))
-end
-
 kron(a::Vector, b::Vector) = [ a[i]*b[j] | i=1:length(a), j=1:length(b) ]
 
 function kron{T,S}(a::Matrix{T}, b::Matrix{S})

j/sort.j

@@ -313,3 +313,41 @@ function issorted(v::AbstractVector)
   end
   return true
 end
+
+function _jl_quickselect(a::AbstractVector, k::Integer, lo::Integer, hi::Integer)
+    if k < lo || k > hi; error("k is out of bounds"); end
+
+    while true
+
+        if lo == hi; return a[lo]; end
+
+        # pivot_ind = partition(a, lo, hi)
+        i, j = lo, hi
+        pivot = a[randival(lo,hi)]
+        while i < j
+            while a[i] < pivot; i += 1; end
+            while a[j] > pivot; j -= 1; end
+            if a[i] == a[j]
+                i += 1
+            elseif i < j
+                a[i], a[j] = a[j], a[i]
+            end
+        end
+        pivot_ind = j
+
+        length = pivot_ind - lo + 1
+        if k == length
+            return a[pivot_ind]
+        elseif k <  length
+            hi = pivot_ind - 1
+        else
+            lo = pivot_ind + 1
+            k = k - length
+        end
+
+    end # while true...
+
+end
+
+select(a::AbstractVector, k::Integer) = _jl_quickselect(copy(a), k, 1, length(a))
+select!(a::AbstractVector, k::Integer) = _jl_quickselect(a, k, 1, length(a))

j/statistics.j

@@ -0,0 +1,13 @@
+mean(V::AbstractVector) = sum(V) / length(V)
+
+function std(V::AbstractVector)
+    n = numel(V)
+    m = mean(V)
+    s = 0.0
+    for i=1:n
+        s += (V[i] - m)^2
+    end
+    return sqrt(s/(n-1))
+end
+
+median(V::AbstractVector) = select(V, div(numel(V),2))

j/sysimg.j

@@ -98,6 +98,7 @@ load("math.j")
 load("math_libm.j")
 load("sort.j")
 load("combinatorics.j")
+load("statistics.j")
 
 # linear algebra
 load("linalg.j")