diff --git a/src/besselj.jl b/src/besselj.jl
index b2959e7..d8a1843 100644
--- a/src/besselj.jl
+++ b/src/besselj.jl
@@ -57,6 +57,13 @@ function _besselj0(x::Float64)
             q2 = (-1/8, 25/384, -1073/5120, 375733/229376)
             p = evalpoly(x2, p2)
             q = evalpoly(x2, q2)
+            if x > 1e15
+                a = SQ2OPI(T) * sqrt(xinv) * p
+                xn = muladd(xinv, q, -PIO4(T))
+                s_x, c_x = sincos(x)
+                s_xn, c_xn = sincos(xn)
+                return a * (c_x * c_xn - s_x * s_xn)
+            end
         end
 
         a = SQ2OPI(T) * sqrt(xinv) * p
@@ -126,6 +133,13 @@ function _besselj1(x::Float64)
             q2 = (3/8, -21/128, 1899/5120, -543483/229376)
             p = evalpoly(x2, p2)
             q = evalpoly(x2, q2)
+            if x > 1e15
+                a = SQ2OPI(T) * sqrt(xinv) * p
+                xn = muladd(xinv, q, -3 * PIO4(T))
+                s_x, c_x = sincos(x)
+                s_xn, c_xn = sincos(xn)
+                return s * a * (c_x * c_xn - s_x * s_xn)
+            end
         end
 
         a = SQ2OPI(T) * sqrt(xinv) * p
diff --git a/src/bessely.jl b/src/bessely.jl
index c392b64..adeefa5 100644
--- a/src/bessely.jl
+++ b/src/bessely.jl
@@ -79,10 +79,17 @@ function _bessely0_compute(x::Float64)
             q2 = (-1/8, 25/384, -1073/5120, 375733/229376)
             p = evalpoly(x2, p2)
             q = evalpoly(x2, q2)
+            if x > 1e15
+                a = SQ2OPI(T) * sqrt(xinv) * p
+                xn = muladd(xinv, q, -PIO4(T))
+                s_x, c_x = sincos(x)
+                s_xn, c_xn = sincos(xn)
+                return a * (s_x * c_xn + c_x * s_xn)
+            end
         end
 
         a = SQ2OPI(T) * sqrt(xinv) * p
-        xn = muladd(xinv, q, - PIO4(T))
+        xn = muladd(xinv, q, -PIO4(T))
 
         # the following computes b = sin(x + xn) more accurately
         # see src/misc.jl
@@ -159,6 +166,13 @@ function _bessely1_compute(x::Float64)
             q2 = (3/8, -21/128, 1899/5120, -543483/229376)
             p = evalpoly(x2, p2)
             q = evalpoly(x2, q2)
+            if x > 1e15
+                a = SQ2OPI(T) * sqrt(xinv) * p
+                xn = muladd(xinv, q, - 3 * PIO4(T))
+                s_x, c_x = sincos(x)
+                s_xn, c_xn = sincos(xn)
+                return a * (s_x * c_xn + c_x * s_xn)
+            end
         end
 
         a = SQ2OPI(T) * sqrt(xinv) * p
diff --git a/test/besselj_test.jl b/test/besselj_test.jl
index 27d91e2..5ec2d69 100644
--- a/test/besselj_test.jl
+++ b/test/besselj_test.jl
@@ -72,6 +72,11 @@ j1_32 = besselj1.(Float32.(x))
 @test besselj1(-80.0f0) ≈ SpecialFunctions.besselj1(-80.0f0)
 @test besselj1(-80.0) ≈ SpecialFunctions.besselj1(-80.0)
 
+# tests for very large inputs
+x = [1e12, 5e12, 1e13, 5e13, 1e14, 5e14, 1e15, 5e15, 1e16, 5e16, 1e17, 5e17, 1e18, 5e18, 1e19, 5e19, 1e20, 1e22, 1e25, 1e30, 1e40]
+@test besselj0.(x) ≈ SpecialFunctions.besselj0.(x)
+@test besselj1.(x) ≈ SpecialFunctions.besselj1.(x)
+
 ## Tests for besselj 
 
 #=
diff --git a/test/bessely_test.jl b/test/bessely_test.jl
index e165885..9f70ae9 100644
--- a/test/bessely_test.jl
+++ b/test/bessely_test.jl
@@ -68,6 +68,11 @@ y1_32 = bessely1.(Float32.(x))
 @test bessely1(Inf32) == zero(Float32)
 @test bessely1(Inf64) == zero(Float64)
 
+# tests for very large inputs
+x = [1e12, 5e12, 1e13, 5e13, 1e14, 5e14, 1e15, 5e15, 1e16, 5e16, 1e17, 5e17, 1e18, 5e18, 1e19, 5e19, 1e20, 1e22, 1e25, 1e30, 1e40]
+@test bessely0.(x) ≈ SpecialFunctions.bessely0.(x)
+@test bessely1.(x) ≈ SpecialFunctions.bessely1.(x)
+
 ## Tests for bessely
 
 ## test all numbers and orders for 0<nu<100