Fix Hurwitz zeta function with q=1

include/xsf/zeta.h

@@ -345,11 +345,24 @@ XSF_HOST_DEVICE inline double riemann_zeta(double x) { return cephes::riemann_ze
 
 XSF_HOST_DEVICE inline float riemann_zeta(float x) { return riemann_zeta(static_cast<double>(x)); }
 
-XSF_HOST_DEVICE inline double zeta(double x, double q) { return cephes::zeta(x, q); }
+XSF_HOST_DEVICE inline double zeta(double x, double q) {
+    if (q == 1.0) {
+        return riemann_zeta(x);
+    }
+    return cephes::zeta(x, q);
+}
 
-XSF_HOST_DEVICE inline float zeta(float x, float q) { return zeta(static_cast<double>(x), static_cast<double>(q)); }
+XSF_HOST_DEVICE inline float zeta(float x, float q) {
+    if (q == 1.0f) {
+        return riemann_zeta(x);
+    }
+    return zeta(static_cast<double>(x), static_cast<double>(q));
+}
 
 XSF_HOST_DEVICE inline std::complex<double> zeta(std::complex<double> z, double q) {
+    if (q == 1.0) {
+        return riemann_zeta(z);
+    }
     if (z.imag() == 0.0) {
         return zeta(z.real(), q);
     }

tests/xsf_tests/test_zeta.cpp

@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+#include <catch2/generators/catch_generators_range.hpp>
+#include <tuple>
+#include <xsf/zeta.h>
+
+TEST_CASE("zeta(x, q=1) matches riemann_zeta for all types", "[zeta][xsf_tests]") {
+    SECTION("double") {
+        double x = GENERATE(range(-10.0, 10.0, 0.1));
+        REQUIRE(xsf::zeta(x, 1.0) == xsf::riemann_zeta(x));
+    }
+
+    SECTION("float") {
+        float x = GENERATE(range(-10.0f, 10.0f, 0.5f));
+        REQUIRE(xsf::zeta(x, 1.0f) == xsf::riemann_zeta(x));
+    }
+
+    SECTION("complex") {
+        using std::complex;
+        double re = GENERATE(range(0.5, 5.0, 0.5));
+        double im = GENERATE(range(-2.0, 2.0, 0.5));
+        complex<double> z(re, im);
+        REQUIRE(xsf::zeta(z, 1.0) == xsf::riemann_zeta(z));
+    }
+}