diff --git a/base/math.jl b/base/math.jl
index 3e87f55168a5e..ba23840267f4c 100644
--- a/base/math.jl
+++ b/base/math.jl
@@ -536,9 +536,6 @@ quadrant of the return value.
 atan2(y::Real, x::Real) = atan2(promote(float(y),float(x))...)
 atan2(y::T, x::T) where {T<:AbstractFloat} = Base.no_op_err("atan2", T)
 
-atan2(y::Float64, x::Float64) = ccall((:atan2,libm), Float64, (Float64, Float64,), y, x)
-atan2(y::Float32, x::Float32) = ccall((:atan2f,libm), Float32, (Float32, Float32), y, x)
-
 max(x::T, y::T) where {T<:AbstractFloat} = ifelse((y > x) | (signbit(y) < signbit(x)),
                                     ifelse(isnan(x), x, y), ifelse(isnan(y), y, x))
 
diff --git a/base/special/rem_pio2.jl b/base/special/rem_pio2.jl
index 9d8f830c03f92..0f45fc0fd0db9 100644
--- a/base/special/rem_pio2.jl
+++ b/base/special/rem_pio2.jl
@@ -58,6 +58,7 @@ Return the high word of `x` as a `UInt32`.
 Return positive part of the high word of `x` as a `UInt32`.
 """
 @inline poshighword(x::UInt64) = unsafe_trunc(UInt32,x >> 32)&0x7fffffff
+@inline poshighword(x::Float32) = reinterpret(UInt32, x)&0x7fffffff
 @inline poshighword(x::Float64) = poshighword(reinterpret(UInt64, x))
 
 @inline function cody_waite_2c_pio2(x::Float64, fn, n)
diff --git a/base/special/trig.jl b/base/special/trig.jl
index a9c7c7a699a46..9da7459d4ea46 100644
--- a/base/special/trig.jl
+++ b/base/special/trig.jl
@@ -175,6 +175,107 @@ function asin(x::T) where T<:Union{Float32, Float64}
     return asin_kernel(t, x)
 end
 
+# atan2 methods
+ATAN2_PI_LO(::Type{Float32}) = -8.7422776573f-08
+ATAN2_RATIO_BIT_SHIFT(::Type{Float32}) = 23
+ATAN2_RATIO_THRESHOLD(::Type{Float32}) = 26
+
+ATAN2_PI_LO(::Type{Float64}) = 1.2246467991473531772E-16
+ATAN2_RATIO_BIT_SHIFT(::Type{Float64}) = 20
+ATAN2_RATIO_THRESHOLD(::Type{Float64}) = 60
+
+function atan2(y::T, x::T) where T<:Union{Float32, Float64}
+    # Method :
+    #    M1) Reduce y to positive by atan2(y,x)=-atan2(-y,x).
+    #    M2) Reduce x to positive by (if x and y are unexceptional):
+    #        ARG (x+iy) = arctan(y/x)          ... if x > 0,
+    #        ARG (x+iy) = pi - arctan[y/(-x)]   ... if x < 0,
+    #
+    # Special cases:
+    #
+    #    S1) ATAN2((anything), NaN ) is NaN;
+    #    S2) ATAN2(NAN , (anything) ) is NaN;
+    #    S3) ATAN2(+-0, +(anything but NaN)) is +-0  ;
+    #    S4) ATAN2(+-0, -(anything but NaN)) is +-pi ;
+    #    S5) ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2;
+    #    S6) ATAN2(+-(anything but INF and NaN), +INF) is +-0 ;
+    #    S7) ATAN2(+-(anything but INF and NaN), -INF) is +-pi;
+    #    S8) ATAN2(+-INF,+INF ) is +-pi/4 ;
+    #    S9) ATAN2(+-INF,-INF ) is +-3pi/4;
+    #    S10) ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2;
+    if isnan(x) || isnan(y) # S1 or S2
+        return T(NaN)
+    end
+
+    if x == T(1.0) # then y/x = y and x > 0, see M2
+        return atan(y)
+    end
+    # generate an m ∈ {0, 1, 2, 3} to branch off of
+    m = 2*signbit(x) + 1*signbit(y)
+
+    if iszero(y)
+        if m == 0 || m == 1
+            return y # atan(+-0, +anything) = +-0
+        elseif m == 2
+            return T(pi) # atan(+0, -anything) = pi
+        elseif m == 3
+            return -T(pi) # atan(-0, -anything) =-pi
+        end
+    elseif iszero(x)
+        return flipsign(T(pi)/2, y)
+    end
+
+    if isinf(x)
+        if isinf(y)
+            if m == 0
+                return T(pi)/4  # atan(+Inf), +Inf))
+            elseif m == 1
+                return -T(pi)/4 # atan(-Inf), +Inf))
+            elseif m == 2
+                return 3*T(pi)/4 # atan(+Inf, -Inf)
+            elseif m == 3
+                return -3*T(pi)/4 # atan(-Inf,-Inf)
+            end
+        else
+            if m == 0
+                return zero(T)  # atan(+...,+Inf) */
+            elseif m == 1
+                return -zero(T) # atan(-...,+Inf) */
+            elseif m == 2
+                return T(pi)    # atan(+...,-Inf) */
+            elseif m == 3
+                return -T(pi)   # atan(-...,-Inf) */
+            end
+        end
+    end
+
+    # x wasn't Inf, but y is
+    isinf(y) && return copysign(T(pi)/2, y)
+
+    ypw = poshighword(y)
+    xpw = poshighword(x)
+    # compute y/x for Float32
+    k = reinterpret(Int32, ypw-xpw)>>ATAN2_RATIO_BIT_SHIFT(T)
+
+    if k > ATAN2_RATIO_THRESHOLD(T) # |y/x| >  threshold
+        z=T(pi)/2+T(0.5)*ATAN2_PI_LO(T)
+        m&=1;
+    elseif x<0 && k < -ATAN2_RATIO_THRESHOLD(T) # 0 > |y|/x > threshold
+        z = zero(T)
+    else #safe to do y/x
+        z = atan(abs(y/x))
+    end
+
+    if m == 0
+        return z # atan(+,+)
+    elseif m == 1
+        return -z # atan(-,+)
+    elseif m == 2
+        return T(pi)-(z-ATAN2_PI_LO(T)) # atan(+,-)
+    else # default case m == 3
+        return (z-ATAN2_PI_LO(T))-T(pi) # atan(-,-)
+    end
+end
 # acos methods
 ACOS_X_MIN_THRESHOLD(::Type{Float32}) = 2.0f0^-26
 ACOS_X_MIN_THRESHOLD(::Type{Float64}) = 2.0^-57
diff --git a/test/math.jl b/test/math.jl
index c0645ba4cfbfc..22cd551fedb25 100644
--- a/test/math.jl
+++ b/test/math.jl
@@ -704,6 +704,60 @@ end
     end
 end
 
+@testset "atan2" begin
+    for T in (Float32, Float64)
+        @test atan2(T(NaN), T(NaN)) === T(NaN)
+        @test atan2(T(NaN), T(0.1)) === T(NaN)
+        @test atan2(T(0.1), T(NaN)) === T(NaN)
+        r = T(randn())
+        absr = abs(r)
+        # y zero
+        @test atan2(T(r), one(T)) === atan(T(r))
+        @test atan2(zero(T), absr) === zero(T)
+        @test atan2(-zero(T), absr) === -zero(T)
+        @test atan2(zero(T), -absr) === T(pi)
+        @test atan2(-zero(T), -absr) === -T(pi)
+        # x zero and y not zero
+        @test atan2(one(T), zero(T)) === T(pi)/2
+        @test atan2(-one(T), zero(T)) === -T(pi)/2
+        # isinf(x) == true && isinf(y) == true
+        @test atan2(T(Inf), T(Inf)) === T(pi)/4 # m == 0 (see atan2 code)
+        @test atan2(-T(Inf), T(Inf)) === -T(pi)/4 # m == 1
+        @test atan2(T(Inf), -T(Inf)) === 3*T(pi)/4 # m == 2
+        @test atan2(-T(Inf), -T(Inf)) === -3*T(pi)/4 # m == 3
+        # isinf(x) == true && isinf(y) == false
+        @test atan2(absr, T(Inf)) === zero(T) # m == 0
+        @test atan2(-absr, T(Inf)) === -zero(T) # m == 1
+        @test atan2(absr, -T(Inf)) === T(pi) # m == 2
+        @test atan2(-absr, -T(Inf)) === -T(pi) # m == 3
+        # isinf(y) == true && isinf(x) == false
+        @test atan2(T(Inf), absr) === T(pi)/2
+        @test atan2(-T(Inf), absr) === -T(pi)/2
+        @test atan2(T(Inf), -absr) === T(pi)/2
+        @test atan2(-T(Inf), -absr) === -T(pi)/2
+        # |y/x| above high threshold
+        atanpi = T(1.5707963267948966)
+        @test atan2(T(2.0^61), T(1.0)) === atanpi # m==0
+        @test atan2(-T(2.0^61), T(1.0)) === -atanpi # m==1
+        @test atan2(T(2.0^61), -T(1.0)) === atanpi # m==2
+        @test atan2(-T(2.0^61), -T(1.0)) === -atanpi # m==3
+        @test atan2(-T(Inf), -absr) === -T(pi)/2
+        # |y|/x between 0 and low threshold
+        @test atan2(T(2.0^-61), -T(1.0)) === T(pi) # m==2
+        @test atan2(-T(2.0^-61), -T(1.0)) === -T(pi) # m==3
+        # y/x is "safe" ("arbitrary values", just need to hit the branch)
+        _ATAN2_PI_LO(::Type{Float32}) = -8.7422776573f-08
+        _ATAN2_PI_LO(::Type{Float64}) = 1.2246467991473531772E-16
+        @test atan2(T(5.0), T(2.5)) === atan(abs(T(5.0)/T(2.5)))
+        @test atan2(-T(5.0), T(2.5)) === -atan(abs(-T(5.0)/T(2.5)))
+        @test atan2(T(5.0), -T(2.5)) === T(pi)-(atan(abs(T(5.0)/-T(2.5)))-_ATAN2_PI_LO(T))
+        @test atan2(-T(5.0), -T(2.5)) === -(T(pi)-atan(abs(-T(5.0)/-T(2.5)))-_ATAN2_PI_LO(T))
+        @test atan2(T(1235.2341234), T(2.5)) === atan(abs(T(1235.2341234)/T(2.5)))
+        @test atan2(-T(1235.2341234), T(2.5)) === -atan(abs(-T(1235.2341234)/T(2.5)))
+        @test atan2(T(1235.2341234), -T(2.5)) === T(pi)-(atan(abs(T(1235.2341234)/-T(2.5)))-_ATAN2_PI_LO(T))
+        @test atan2(-T(1235.2341234), -T(2.5)) === -(T(pi)-(atan(abs(-T(1235.2341234)/T(2.5)))-_ATAN2_PI_LO(T)))
+    end
+end
 @testset "acos #23283" begin
     for T in (Float32, Float64)
         @test acos(zero(T)) === T(pi)/2