diff --git a/include/xsf/amos/amos.h b/include/xsf/amos/amos.h
index 815a99b81f..87c9eb36d6 100644
--- a/include/xsf/amos/amos.h
+++ b/include/xsf/amos/amos.h
@@ -4132,7 +4132,7 @@ namespace amos {
                 continue;
             }
             cs = -zr + std::log(s1);
-            cs = (exp(std::real(cs)) / tol) * (cos(std::imag(cs)) + sin(std::imag(cs) * std::complex<double>(0, 1)));
+            cs = (exp(std::real(cs)) / tol) * std::complex<double>(cos(std::imag(cs)), sin(std::imag(cs)));
             if (!uchk(cs, *ascle, tol)) {
                 y[i - 1] = cs;
                 nz -= 1;
diff --git a/tests/xsf_tests/README.md b/tests/xsf_tests/README.md
new file mode 100644
index 0000000000..b85658a899
--- /dev/null
+++ b/tests/xsf_tests/README.md
@@ -0,0 +1 @@
+These tests are independent of SciPy.
diff --git a/tests/xsf_tests/test_cyl_bessel_k.cpp b/tests/xsf_tests/test_cyl_bessel_k.cpp
new file mode 100644
index 0000000000..d1691fbcc0
--- /dev/null
+++ b/tests/xsf_tests/test_cyl_bessel_k.cpp
@@ -0,0 +1,20 @@
+#include "../testing_utils.h"
+#include <complex>
+#include <tuple>
+#include <xsf/bessel.h>
+
+TEST_CASE("cyl_bessel_k gh-46", "[cyl_bessel_k][xsf_tests]") {
+    using test_case = std::tuple<double, std::complex<double>, std::complex<double>, double>;
+    using std::complex;
+    // Reference values were computed with the Python library mpmath.
+    auto [v, z, ref, rtol] = GENERATE(
+        test_case{0.0, complex{680.0, -1000.0}, complex{1.901684871999608e-298, 1.713412341479591e-297}, 1e-13},
+        test_case{0.0, complex{680.0, -680.0}, complex{-4.553730032944803e-298, 1.878727010109855e-297}, 1e-13},
+        test_case{0.0, complex{25.0, 100.0}, complex{1.699267100365868e-12, -2.234890030902166e-13}, 5e-16}
+    );
+    const complex w = xsf::cyl_bessel_k(v, z);
+    const auto rel_error = xsf::extended_relative_error(w, ref);
+
+    CAPTURE(v, z, w, ref, rtol, rel_error);
+    REQUIRE(rel_error <= rtol);
+}