Fix Hurwitz zeta function when calculated using a Bernoulli polynomial and a is nonpositive integer, fixes #50

Author: mtommila

apfloat/src/main/java/org/apfloat/HurwitzZetaHelper.java

@@ -60,6 +60,12 @@ public static Apcomplex zeta(Apcomplex s, Apcomplex a)
              precision = ApfloatHelper.extendPrecision(targetPrecision);
         s = ApfloatHelper.ensurePrecision(s, precision);
         a = ApfloatHelper.ensurePrecision(a, precision);
+        if (isNonPositiveInteger(s))
+        {
+            // Handle cases where the result might be genuinely zero
+            long n1 = ApfloatHelper.longValueExact(one.subtract(s.real().truncate()));
+            return bernoulliB(n1, a).divide(new Apfloat(n1, targetPrecision, radix)).negate();
+        }
         if (isNonPositiveInteger(a))
         {
             if (s.real().signum() < 0 || s.isZero())
@@ -86,12 +92,6 @@ public static Apcomplex zeta(Apcomplex s, Apcomplex a)
 
             throw new ArithmeticException("Zeta of second argument nonpositive integer");
         }
-        if (isNonPositiveInteger(s))
-        {
-            // Handle cases where the result might be genuinely zero
-            long n1 = ApfloatHelper.longValueExact(one.subtract(s.real().truncate()));
-            return bernoulliB(n1, a).divide(new Apint(n1, radix)).negate();
-        }
         if (precision == Apfloat.INFINITE)
         {
             throw new InfiniteExpansionException("Cannot calculate zeta function to infinite precision");
@@ -273,8 +273,13 @@ private static Apfloat expPrecision(Apfloat x)
     // z assumed to be of extended precision already
     private static Apcomplex bernoulliB(long n, Apcomplex z)
     {
+        assert (n > 0);
         long precision = z.precision();
         int radix = z.radix();
+        if (z.isZero())
+        {
+            return AprationalMath.bernoulli(n, radix);
+        }
         Apcomplex sum = Apcomplex.ZEROS[radix],
                   numerator = new Apfloat(1, precision, radix),
                   denominator = numerator;

apfloat/src/test/java/org/apfloat/ApcomplexMathTest.java

@@ -3297,6 +3297,14 @@ public static void testZetaHurwitz()
         assertEquals("-1+2i,-3 precision", 6, a.precision());
         assertEquals("-1+2i,-3 value", new Apcomplex("(210.002,2354.18)"), a, new Apfloat("5e-2"));
 
+        a = ApcomplexMath.zeta(new Apcomplex("0"), new Apcomplex("-1.00000"));
+        assertEquals("0,-1 precision", 6, a.precision());
+        assertEquals("0,-1 value", new Apcomplex("(1.50000)"), a, new Apfloat("5e-5"));
+
+        a = ApcomplexMath.zeta(new Apcomplex("0"), new Apcomplex("-2.00000"));
+        assertEquals("0,-2 precision", 6, a.precision());
+        assertEquals("0,-2 value", new Apcomplex("(2.50000)"), a, new Apfloat("5e-5"));
+
         a = ApcomplexMath.zeta(new Apcomplex("(2,3)"), new Apcomplex("(4,5)"));
         assertEquals("2+3i,4+5i precision", 1, a.precision());
         assertEquals("2+3i,4+5i value", new Apcomplex("(0.2,-0.9)"), a, new Apfloat("5e-1"));
@@ -5939,6 +5947,10 @@ public static void testBernoulliB()
         assertEquals("4, -3 precision", 6, a.precision());
         assertEquals("4, -3 value", new Apcomplex("143.967"), a, new Apfloat("5e-3"));
 
+        a = ApcomplexMath.bernoulliB(1, new Apcomplex("-2.00000"));
+        assertEquals("1, -2 precision", 6, a.precision());
+        assertEquals("1, -2 value", new Apcomplex("-2.50000"), a, new Apfloat("5e-5"));
+
         a = ApcomplexMath.bernoulliB(3, new Apcomplex(new Apfloat("0.1", 18, 2), new Apfloat("0.11", 18, 2)));
         assertEquals("3, 0.5 + 0.75i precision", 18, a.precision());
         assertEquals("3, 0.5 + 0.75i radix", 2, a.radix());

apfloat/src/test/java/org/apfloat/ApfloatMathTest.java

@@ -3125,7 +3125,7 @@ public static void testZetaHurwitz()
         assertEquals("-3,-3 value", new Apfloat("-35.9917"), a, new Apfloat("5e-4"));
 
         a = ApfloatMath.zeta(new Apfloat("-3.00000"), new Apfloat("-1.00000"));
-        assertEquals("-3,-1 precision", 5, a.precision());
+        assertEquals("-3,-1 precision", 6, a.precision());
         assertEquals("-3,-1 value", new Apfloat("-0.991667"), a, new Apfloat("5e-6"));
 
         a = ApfloatMath.zeta(new Apfloat("-3.00000"), new Apfloat("0"));