Some fixes to overflow handling

apfloat/src/main/java/org/apfloat/ApcomplexMath.java

@@ -1737,15 +1737,7 @@ public static Apcomplex gamma(Apcomplex z)
                 {
                     throw new ArithmeticException("Gamma of negative integer");
                 }
-                long n;
-                try
-                {
-                    n = z.real().longValueExact();
-                }
-                catch (ArithmeticException ae)
-                {
-                    throw new OverflowException("Overflow");
-                }
+                long n = ApfloatHelper.longValueExact(z.real().truncate());
                 // If n is extremely large and precision is relatively low, then computing it as gamma is faster
                 // If n is relatively small compared to precision, then computing it as factorial is faster
                 double gammaEffort = Math.log(precision) * precision * precision;
@@ -1987,7 +1979,7 @@ public static Apcomplex logGamma(Apcomplex z)
         Iterator<Aprational> bernoulli2 = AprationalMath.bernoullis2(n, radix);
         for (long k = 1; k <= n; k++)
         {
-            long k2 = Math.multiplyExact(k,  2);
+            long k2 = Util.multiplyExact(k,  2);
             Apcomplex term = bernoulli2.next().precision(workingPrecision).divide(new Apint(k2, radix).multiply(new Apint(k2 - 1, radix)).multiply(zp));
             if (k < n)
             {
@@ -2220,7 +2212,7 @@ public static Apcomplex digamma(Apcomplex z)
         Iterator<Aprational> bernoulli2 = AprationalMath.bernoullis2(n, radix);
         for (long k = 1; k <= n; k++)
         {
-            long k2 = Math.multiplyExact(k,  2);
+            long k2 = Util.multiplyExact(k,  2);
             zp = zp.multiply(z2);
             Apcomplex term = bernoulli2.next().precision(precision).divide(new Apint(k2, radix).multiply(zp));
 

apfloat/src/main/java/org/apfloat/ApfloatHelper.java

@@ -1,7 +1,7 @@
 /*
  * MIT License
  *
- * Copyright (c) 2002-2023 Mikko Tommila
+ * Copyright (c) 2002-2024 Mikko Tommila
  *
  * Permission is hereby granted, free of charge, to any person obtaining a copy
  * of this software and associated documentation files (the "Software"), to deal
@@ -41,7 +41,7 @@
 /**
  * Various utility methods related to apfloats.
  *
- * @version 1.10.1
+ * @version 1.13.0
  * @author Mikko Tommila
  */
 
@@ -801,6 +801,20 @@ public int read(byte[] b, int off, int len)
         return new PushbackReader(new InputStreamReader(in, "ISO-8859-1"));
     }
 
+    // Converts ArithmeticException to OverflowException if the Apint overflows a long
+    public static long longValueExact(Apint x)
+        throws OverflowException
+    {
+        try
+        {
+            return x.longValueExact();
+        }
+        catch (ArithmeticException ae)
+        {
+            throw new OverflowException("Overflow", ae);
+        }
+    }
+
     private static int getDefaultRadix()
         throws NumberFormatException
     {

apfloat/src/main/java/org/apfloat/ApfloatMath.java

@@ -2420,7 +2420,7 @@ else if (precision == Apfloat.INFINITE)
                 four = new Apfloat(4, workingPrecision, radix),
                 log4 = log(four),
                 pi = pi(workingPrecision, radix),
-                d = pow(two, Math.multiplyExact(2,  n) - 1),
+                d = pow(two, Util.multiplyExact(2,  n) - 1),
                 dnk = d,
                 z = Apfloat.ZERO;
         for (long k = n - 1; k >= 0; k--)
@@ -2502,7 +2502,7 @@ else if (precision == Apfloat.INFINITE)
         Iterator<Aprational> bernoullis2 = AprationalMath.bernoullis2Small(radix);
         for (long n = 1; ; n++)
         {
-            f = f.multiply(twopi2).divide(new Apint(Math.multiplyExact(2, n) - 1, radix).multiply(new Apint(2 * n, radix)));
+            f = f.multiply(twopi2).divide(new Apint(Util.multiplyExact(2, n) - 1, radix).multiply(new Apint(2 * n, radix)));
             Apfloat z = abs(bernoullis2.next()).multiply(f).subtract(one);
             if (z.scale() < -precision)
             {

apfloat/src/main/java/org/apfloat/ApintMath.java

@@ -1,7 +1,7 @@
 /*
  * MIT License
  *
- * Copyright (c) 2002-2023 Mikko Tommila
+ * Copyright (c) 2002-2024 Mikko Tommila
  *
  * Permission is hereby granted, free of charge, to any person obtaining a copy
  * of this software and associated documentation files (the "Software"), to deal
@@ -31,7 +31,7 @@
 /**
  * Various mathematical functions for arbitrary precision integers.
  *
- * @version 1.11.0
+ * @version 1.13.0
  * @author Mikko Tommila
  */
 
@@ -708,7 +708,14 @@ public static Apint binomial(long n, long k, int radix)
         {
             if (k >= 0)
             {
-                n = Math.subtractExact(k, n) - 1;
+                try
+                {
+                    n = Math.subtractExact(k, n) - 1;
+                }
+                catch (ArithmeticException ae)
+                {
+                    return binomial(new Apint(n, radix), new Apint(k, radix));
+                }
             }
             else if (k <= n)
             {
@@ -724,7 +731,14 @@ else if (k <= n)
         }
         else if (k < 0)
         {
-            k = Math.subtractExact(n, k);
+            try
+            {
+                k = Math.subtractExact(n, k);
+            }
+            catch (ArithmeticException ae)
+            {
+                return binomial(new Apint(n, radix), new Apint(k, radix));
+            }
         }
         if (k < 0 || k > n)
         {
@@ -802,7 +816,8 @@ else if (k.signum() < 0)
             // Optimize performance
             k = n.subtract(k);
         }
-        Apint b = pochhammer(n.subtract(k).add(one), k).divide(factorial(k.longValueExact(), radix));
+        Apint f = factorial(ApfloatHelper.longValueExact(k), radix),
+              b = pochhammer(n.subtract(k).add(one), k).divide(f);
         return (negate ? b.negate() : b);
     }
 

apfloat/src/main/java/org/apfloat/AprationalMath.java

@@ -1,7 +1,7 @@
 /*
  * MIT License
  *
- * Copyright (c) 2002-2023 Mikko Tommila
+ * Copyright (c) 2002-2024 Mikko Tommila
  *
  * Permission is hereby granted, free of charge, to any person obtaining a copy
  * of this software and associated documentation files (the "Software"), to deal
@@ -500,8 +500,9 @@ public static Aprational binomial(Aprational n, Aprational k)
         {
             return Apint.ZEROS[radix];
         }
-        Apint one = new Apint(1, radix);
-        return pochhammer(n.subtract(k).add(one), k.numerator()).divide(ApintMath.factorial(k.longValueExact(), radix));
+        Apint one = new Apint(1, radix),
+              f = ApintMath.factorial(ApfloatHelper.longValueExact(k.truncate()), radix);
+        return pochhammer(n.subtract(k).add(one), k.numerator()).divide(f);
     }
 
     /**
@@ -635,7 +636,7 @@ static Aprational bernoulliBig(long n, int radix)
         Apint one = Apint.ONES[radix],
               two = new Apint(2, radix);
         long p = Math.max(1, (long) Math.ceil(n * Math.log(n) / Math.log(radix))),
-             precision = ApfloatHelper.extendPrecision(Math.multiplyExact(Math.multiplyExact(2, n) + 1, p));  // 2 * n * p + 2 is not enough!
+             precision = ApfloatHelper.extendPrecision(Util.multiplyExact(Util.multiplyExact(2, n) + 1, p));  // 2 * n * p + 2 is not enough!
         Apint f2n1 = ApintMath.factorial(2 * n - 1,  radix);
         Apfloat z = ApfloatMath.scale(new Apfloat(1, precision, radix), -p),
                 v = ApfloatMath.scale(f2n1.multiply(ApfloatMath.tan(z)), -p);
@@ -781,7 +782,7 @@ static Iterator<Aprational> bernoullis2Big(long n, int radix)
         Apint one = Apint.ONES[radix],
               two = new Apint(2, radix);
         long p = (long) Math.ceil(n * Math.log(n) / Math.log(radix)),
-             precision = ApfloatHelper.extendPrecision(Math.multiplyExact(Math.multiplyExact(2, n) + 1, p));  // 2 * n * p + 2 is not enough!
+             precision = ApfloatHelper.extendPrecision(Util.multiplyExact(Util.multiplyExact(2, n) + 1, p));  // 2 * n * p + 2 is not enough!
         Apint f2n1 = ApintMath.factorial(2 * n - 1,  radix);
         Apfloat z = ApfloatMath.scale(new Apfloat(1, precision, radix), -p);
 

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

@@ -1,7 +1,7 @@
 /*
  * MIT License
  *
- * Copyright (c) 2002-2023 Mikko Tommila
+ * Copyright (c) 2002-2024 Mikko Tommila
  *
  * Permission is hereby granted, free of charge, to any person obtaining a copy
  * of this software and associated documentation files (the "Software"), to deal
@@ -60,7 +60,7 @@ public static Apcomplex zeta(Apcomplex s, Apcomplex a)
             {
                 // Use recurrence formula: zeta(s, a) = a^-s + zeta(s, a + 1)
                 Apcomplex t = Apcomplex.ZERO;
-                long i = a.longValueExact();
+                long i = ApfloatHelper.longValueExact(a.real().truncate());
                 while (i++ <= 0)
                 {
                     t = t.add(ApcomplexMath.pow(a, s.negate()));
@@ -110,7 +110,7 @@ private static Apint binarySearch(Apcomplex s, Apcomplex a, long precision, long
               two = new Apint(2, radix),
               min = ApintMath.max(two.subtract(a.real().truncate()), one.subtract(s.real()).divide(two).truncate().add(one)),
               N = ApintMath.max(new Apint(precision, radix), min);
-        long M = N.longValueExact(),
+        long M = ApfloatHelper.longValueExact(N),
              low = M,
              high = M;
         Apfloat R = R(s, a, N, M);
@@ -120,7 +120,7 @@ private static Apint binarySearch(Apcomplex s, Apcomplex a, long precision, long
         {
             low = M;
             // R is too large, M is too small
-            M = Math.multiplyExact(2, M);
+            M = Util.multiplyExact(2, M);
             N = new Apint(M, radix);
             R = R(s, a, N, M);
             high = M;

apfloat/src/main/java/org/apfloat/IncompleteGammaHelper.java

@@ -1,7 +1,7 @@
 /*
  * MIT License
  *
- * Copyright (c) 2002-2023 Mikko Tommila
+ * Copyright (c) 2002-2024 Mikko Tommila
  *
  * Permission is hereby granted, free of charge, to any person obtaining a copy
  * of this software and associated documentation files (the "Software"), to deal
@@ -27,11 +27,13 @@
 import java.util.function.BiFunction;
 import java.util.function.LongFunction;
 
+import org.apfloat.spi.Util;
+
 /**
  * Helper class for the incomplete gamma function.
  *
  * @since 1.10.0
- * @version 1.11.1
+ * @version 1.13.0
  * @author Mikko Tommila
  */
 
@@ -103,33 +105,33 @@ private static enum ContinuedFraction
             @Override
             protected long doGetMinIterations(Apcomplex a, Apcomplex z)
             {
-                return (a.real().signum() >= 0 ? 0 : Math.subtractExact(4, Math.multiplyExact(2, a.real().truncate().longValueExact())));
+                return (a.real().signum() >= 0 ? 0 : Util.subtractExact(4, Util.multiplyExact(2, ApfloatHelper.longValueExact(a.real().truncate()))));
             }
         },
         LOWER2(ContinuedFractionType.LOWER, IncompleteGammaHelper::lowerGammaSequenceAlternative)
         {
             @Override
             protected long doGetMinIterations(Apcomplex a, Apcomplex z)
             {
-                return Math.max(a.real().signum() >= 0 ? 0 : Math.subtractExact(3, a.real().truncate().longValueExact()),
-                                Math.subtractExact(2, Math.addExact(a.real().truncate().longValueExact(), z.real().truncate().longValueExact())));
+                return Math.max(a.real().signum() >= 0 ? 0 : Util.subtractExact(3, ApfloatHelper.longValueExact(a.real().truncate())),
+                                Util.subtractExact(2, Util.addExact(a.real().truncate().longValueExact(), ApfloatHelper.longValueExact(z.real().truncate()))));
             }
         },
         UPPER1(ContinuedFractionType.UPPER, IncompleteGammaHelper::upperGammaSequence)
         {
             @Override
             protected long doGetMinIterations(Apcomplex a, Apcomplex z)
             {
-                return Math.max(a.real().signum() <= 0 ? 0 : Math.addExact(2, a.real().truncate().longValueExact()),
-                                Math.addExact(1, Math.subtractExact(a.real().truncate().longValueExact(), z.real().truncate().longValueExact()) / 2));
+                return Math.max(a.real().signum() <= 0 ? 0 : Util.addExact(2, ApfloatHelper.longValueExact(a.real().truncate())),
+                                Util.addExact(1, Util.subtractExact(ApfloatHelper.longValueExact(a.real().truncate()), ApfloatHelper.longValueExact(z.real().truncate())) / 2));
             }
         },
         UPPER2(ContinuedFractionType.UPPER, IncompleteGammaHelper::upperGammaSequenceAlternative)
         {
             @Override
             protected long doGetMinIterations(Apcomplex a, Apcomplex z)
             {
-                return (a.real().signum() <= 0 ? 0 : Math.addExact(Math.multiplyExact(2, a.real().truncate().longValueExact()), 2));
+                return (a.real().signum() <= 0 ? 0 : Util.addExact(Util.multiplyExact(2, ApfloatHelper.longValueExact(a.real().truncate())), 2));
             }
         };
 
@@ -289,15 +291,7 @@ private static GammaValue upperGamma(Apcomplex a, Apcomplex z)
                 // Note that this transformation may be extremely slow if n is large
                 // gamma(-n,z) = (-1)^n/n! gamma(0,z) - e^-z sum[z^(k-n-1)/(-n)_k,{k,1,n}]
                 // See https://functions.wolfram.com/GammaBetaErf/Gamma2/17/02/01/
-                long n;
-                try
-                {
-                    n = a.real().longValueExact(); // If this overflows then the value would overflow anyways
-                }
-                catch (ArithmeticException ae)
-                {
-                    throw new OverflowException(ae.getMessage(), ae);
-                }
+                long n = ApfloatHelper.longValueExact(a.real().truncate()); // If this overflows then the value would overflow anyways
                 return new GammaValue(a, upperGamma(n, z), false);
             }
             algorithms = ContinuedFraction.upperValues();
@@ -591,7 +585,7 @@ private static ContinuedFractionResult continuedFraction(Sequence s, int radix,
         long maxPrecision = 0;
         do
         {
-            n = Math.addExact(n, 1);
+            n = Util.addExact(n, 1);
             an = s.a(n).precision(workingPrecision);
             bn = s.b(n).precision(workingPrecision);
             d = d.multiply(an).add(bn);

apfloat/src/main/java/org/apfloat/ZetaHelper.java

@@ -1,7 +1,7 @@
 /*
  * MIT License
  *
- * Copyright (c) 2002-2023 Mikko Tommila
+ * Copyright (c) 2002-2024 Mikko Tommila
  *
  * Permission is hereby granted, free of charge, to any person obtaining a copy
  * of this software and associated documentation files (the "Software"), to deal
@@ -23,11 +23,13 @@
  */
 package org.apfloat;
 
+import org.apfloat.spi.Util;
+
 /**
  * Helper class for the Riemann zeta function.
  *
  * @since 1.11.0
- * @version 1.11.0
+ * @version 1.13.0
  * @author Mikko Tommila
  */
 
@@ -108,7 +110,7 @@ public static Apcomplex alternatingSum(Apcomplex s)
         double t = Math.abs(s.imag().doubleValue()),
                s12 = ApcomplexMath.abs(one.subtract(ApcomplexMath.pow(two, one.subtract(s.precision(doublePrecision))))).doubleValue();
         long n = (long) ((workingPrecision * Math.log(radix) + t * Math.PI / 2 + Math.log((1 + 2 * t) / s12)) / Math.log(3 + Math.sqrt(8))),
-             n2 = Math.multiplyExact(n, 2);
+             n2 = Util.multiplyExact(n, 2);
         Apfloat denominator = ApfloatMath.factorial(n2, workingPrecision, radix),
                 numerator = new Apint(n, radix).multiply(denominator).multiply(ApfloatMath.pow(four.precision(workingPrecision), n)),
                 d = Apfloat.ZERO;