diff --git a/base/math.jl b/base/math.jl
index 7b1da63c3f625..9df60ed828544 100644
--- a/base/math.jl
+++ b/base/math.jl
@@ -166,15 +166,19 @@ for f in (:atan2, :hypot)
     end
 end
 
-max{T<:FloatingPoint}(x::T, y::T) = ifelse((y > x) | (x != x) | (signbit(y) < signbit(x)), y, x)
+max{T<:FloatingPoint}(x::T, y::T) = ifelse((y > x) | (signbit(y) < signbit(x)),
+                                    ifelse(isnan(y), x, y), ifelse(isnan(x), y, x))
+
 @vectorize_2arg Real max
 
-min{T<:FloatingPoint}(x::T, y::T) = ifelse((y < x) | (x != x) | (signbit(y) > signbit(x)), y, x)
+min{T<:FloatingPoint}(x::T, y::T) = ifelse((y < x) | (signbit(y) > signbit(x)),
+                                    ifelse(isnan(y), x, y), ifelse(isnan(x), y, x))
 @vectorize_2arg Real min
 
-minmax{T<:FloatingPoint}(x::T, y::T) =  ifelse((y < x) | (signbit(y) > signbit(x)), (y, x),
-                                        ifelse((y > x) | (signbit(y) < signbit(x)), (x, y),
-                                        ifelse(x == x, (x, x), (y, y))))
+minmax{T<:FloatingPoint}(x::T, y::T) = ifelse(isnan(x-y), ifelse(isnan(x), (y, y), (x, x)),
+                                       ifelse((y < x) | (signbit(y) > signbit(x)), (y, x),
+                                       ifelse((y > x) | (signbit(y) < signbit(x)), (x, y),
+                                       ifelse(x == x, (x, x), (y, y)))))
 
 ldexp(x::Float64,e::Integer) = ccall((:scalbn,libm),  Float64, (Float64,Int32), x, Int32(e))
 ldexp(x::Float32,e::Integer) = ccall((:scalbnf,libm), Float32, (Float32,Int32), x, Int32(e))
diff --git a/test/numbers.jl b/test/numbers.jl
index b9892923a46bc..6b3e94b33be9c 100644
--- a/test/numbers.jl
+++ b/test/numbers.jl
@@ -52,6 +52,7 @@
 @test 2.0 * 3.0 == 6.
 @test min(1.0,1) == 1
 
+# min, max and minmax
 @test minmax(5, 3) == (3, 5)
 @test minmax(3., 5.) == (3., 5.)
 @test minmax(5., 3.) == (3., 5.)
@@ -61,6 +62,22 @@
 @test min(-0.0,0.0) === min(0.0,-0.0)
 @test max(-0.0,0.0) === max(0.0,-0.0)
 @test minmax(-0.0,0.0) === minmax(0.0,-0.0)
+@test max(-3.2, 5.1) == max(5.1, -3.2) == 5.1
+@test min(-3.2, 5.1) == min(5.1, -3.2) == -3.2
+@test max(-3.2, Inf) == max(Inf, -3.2) == Inf
+@test max(-3.2, NaN) == max(NaN, -3.2) == -3.2
+@test min(5.1, Inf) == min(Inf, 5.1) == 5.1
+@test min(5.1, -Inf) == min(-Inf, 5.1) == -Inf
+@test min(5.1, NaN) == min(NaN, 5.1) == 5.1
+@test min(5.1, -NaN) == min(-NaN, 5.1) == 5.1
+@test minmax(-3.2, 5.1) == (min(-3.2, 5.1), max(-3.2, 5.1))
+@test minmax(-3.2, Inf) == (min(-3.2, Inf), max(-3.2, Inf))
+@test minmax(-3.2, NaN) == (min(-3.2, NaN), max(-3.2, NaN))
+@test (max(Inf,NaN), max(-Inf,NaN), max(Inf,-NaN), max(-Inf,-NaN)) == (Inf, -Inf, Inf, -Inf)
+@test (max(NaN,Inf), max(NaN,-Inf), max(-NaN,Inf), max(-NaN,-Inf)) == (Inf, -Inf, Inf, -Inf)
+@test (min(Inf,NaN), min(-Inf,NaN), min(Inf,-NaN), min(-Inf,-NaN)) == (Inf, -Inf, Inf, -Inf)
+@test (min(NaN,Inf), min(NaN,-Inf), min(-NaN,Inf), min(-NaN,-Inf)) == (Inf, -Inf, Inf, -Inf)
+@test minmax(-Inf,NaN) == (min(-Inf,NaN), max(-Inf,NaN))
 
 # fma
 let x = Int64(7)^7