This is an automated email from the ASF dual-hosted git repository.
asf-gitbox-commits pushed a commit to branch master
in repository https://gitbox.apache.org/repos/asf/commons-numbers.git
The following commit(s) were added to refs/heads/master by this push:
new 8ac25c90 NUMBERS-209: Compute Trigamma on negative arguments using
reflection formula
8ac25c90 is described below
commit 8ac25c90936d47a90b41126dffd0a0511094acf2
Author: Alex Herbert <[email protected]>
AuthorDate: Thu Aug 20 16:35:20 2026 +0100
NUMBERS-209: Compute Trigamma on negative arguments using reflection
formula
---
.../org/apache/commons/numbers/gamma/Trigamma.java | 17 ++++
.../apache/commons/numbers/gamma/TrigammaTest.java | 112 ++++++++++++---------
src/changes/changes.xml | 5 +
3 files changed, 86 insertions(+), 48 deletions(-)
diff --git
a/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/Trigamma.java
b/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/Trigamma.java
index 51fe07ab..be9cd1df 100644
---
a/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/Trigamma.java
+++
b/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/Trigamma.java
@@ -51,6 +51,23 @@ public final class Trigamma {
return x;
}
+ if (x < 0) {
+ // Use the reflection formula:
+ // trigamma(x) + trigamma(1 - x) = pi^2 / sin^2(pi * x)
+ // to fall back into positive values. Without this the increment
+ // loop below does not terminate for large negative x: when
+ // ulp(x) >= 2 the update x += 1 leaves x unchanged (or advances
+ // a single step and then sticks) and x < C_LIMIT holds forever.
+
+ // negative integers are poles
+ if (Math.rint(x) == x) {
+ return Double.POSITIVE_INFINITY;
+ }
+
+ final double s = Math.sin(Math.PI * x);
+ return Math.PI * Math.PI / (s * s) - value(1 - x);
+ }
+
if (x > 0 && x <= S_LIMIT) {
return 1 / (x * x);
}
diff --git
a/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/TrigammaTest.java
b/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/TrigammaTest.java
index 09b7c2d5..be916fc6 100644
---
a/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/TrigammaTest.java
+++
b/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/TrigammaTest.java
@@ -18,62 +18,78 @@ package org.apache.commons.numbers.gamma;
import org.junit.jupiter.api.Assertions;
import org.junit.jupiter.api.Test;
+import org.junit.jupiter.params.ParameterizedTest;
+import org.junit.jupiter.params.provider.CsvSource;
/**
* Tests for {@link Trigamma}.
*/
class TrigammaTest {
- @Test
- void testTrigamma() {
- final double eps = 1e-9; // Allowed relative error.
+ @ParameterizedTest
+ @CsvSource({
+ // Negative integers are poles of the trigamma function
+ "-1e300, Infinity",
+ "-1e18, Infinity",
+ "-1e16, Infinity",
+ "-0x1p53, Infinity",
+ "-63, Infinity",
+ "-62, Infinity",
+ "-2, Infinity",
+ "-1, Infinity",
// computed using webMathematica. For example, to compute
trigamma($i) = Polygamma(1, $i), use
//
//
http://functions.wolfram.com/webMathematica/Evaluated.jsp?name=PolyGamma2&plottype=0&vars={%221%22,%22$i%22}&digits=20
- final double[] data = {
- -98765.4321, 10.332673372988805424,
- -100.5, 9.8597034918700861520,
- -50.5, 9.8499971860824842274,
- -20.5, 9.8219943446498794821,
- -10.5, 9.7787577398148123845,
- -5.5, 9.7033198653394003812,
- -2.5, 9.5392466449891237539,
- -0.5, 8.9348022005446793094,
- -1e-1, 101.92253995947720352,
- -1e-2, 10001.669304101071825,
- -1e-3, 1.0000016473414317771e6,
- -1e-4, 1.0000000164517451070e8,
- -1e-5, 1.0000000001644958108e10,
- 1e-11, 1e22,
- 1e-10, 1e20,
- 1e-9, 1.0000000000000000016e18,
- 1e-8, 1.0000000000000001645e16,
- 1e-7, 1.0000000000000164493e14,
- 1e-6, 1.0000000000016449317e12,
- 1e-5, 1.0000000001644910026e10,
- 1e-4, 1.0000000164469368793e8,
- 1e-3, 1.0000016425331958690e6,
- 1e-2, 10001.621213528313220,
- 1e-1, 101.43329915079275882,
- 1, 1.6449340668482264365,
- 1.5, 0.93480220054467930942,
- 2, 0.64493406684822643647,
- 2.5, 0.49035775610023486497,
- 3, 0.39493406684822643647,
- 3.5, 0.33035775610023486497,
- 4, 0.28382295573711532536,
- 4.5, 0.24872510303901037518,
- 5, 0.22132295573711532536,
- 7.5, 0.14261589669670379977,
- 10, 0.10516633568168574612,
- 20, 0.051270822935203119832,
- 50, 0.020201333226697125806,
- 100, 0.010050166663333571395,
- 12345.6789, 0.000081003281325733214110
- };
- for (int i = data.length - 2; i >= 0; i -= 2) {
- final double value = data[i];
- final double expected = data[i + 1];
- Assertions.assertEquals(1, Trigamma.value(value) / expected, eps,
() -> "trigamma " + value);
+ "-98765.4321, 10.332673372988805424",
+ "-100.5, 9.8597034918700861520",
+ "-63.0001, 1.0000000327412051737e8",
+ "-62.0001, 1.0000000327386856554e8",
+ "-50.5, 9.8499971860824842274",
+ "-20.5, 9.8219943446498794821",
+ "-10.5, 9.7787577398148123845",
+ "-5.5, 9.7033198653394003812",
+ "-2.5, 9.5392466449891237539",
+ "-2.0001, 1.0000000289494954257e8",
+ "-1.0001, 1.0000000264497454070e8",
+ "-0.5, 8.9348022005446793094",
+ "-1e-1, 101.92253995947720352",
+ "-1e-2, 10001.669304101071825",
+ "-1e-3, 1.0000016473414317771e6",
+ "-1e-4, 1.0000000164517451070e8",
+ "-1e-5, 1.0000000001644958108e10",
+ "1e-11, 1e22",
+ "1e-10, 1e20",
+ "1e-9, 1.0000000000000000016e18",
+ "1e-8, 1.0000000000000001645e16",
+ "1e-7, 1.0000000000000164493e14",
+ "1e-6, 1.0000000000016449317e12",
+ "1e-5, 1.0000000001644910026e10",
+ "1e-4, 1.0000000164469368793e8",
+ "1e-3, 1.0000016425331958690e6",
+ "1e-2, 10001.621213528313220",
+ "1e-1, 101.43329915079275882",
+ "1, 1.6449340668482264365",
+ "1.5, 0.93480220054467930942",
+ "2, 0.64493406684822643647",
+ "2.5, 0.49035775610023486497",
+ "3, 0.39493406684822643647",
+ "3.5, 0.33035775610023486497",
+ "4, 0.28382295573711532536",
+ "4.5, 0.24872510303901037518",
+ "5, 0.22132295573711532536",
+ "7.5, 0.14261589669670379977",
+ "10, 0.10516633568168574612",
+ "20, 0.051270822935203119832",
+ "50, 0.020201333226697125806",
+ "100, 0.010050166663333571395",
+ "12345.6789, 0.000081003281325733214110",
+ })
+ void testTrigamma(double value, double expected) {
+ final double eps = 1e-9; // Allowed relative error.
+ final double actual = Trigamma.value(value);
+ if (Double.isFinite(expected)) {
+ Assertions.assertEquals(expected, actual, expected * eps, () ->
"trigamma: " + value);
+ } else {
+ Assertions.assertEquals(expected, actual, () -> "trigamma: " +
value);
}
}
diff --git a/src/changes/changes.xml b/src/changes/changes.xml
index 3b60e51a..04a41560 100644
--- a/src/changes/changes.xml
+++ b/src/changes/changes.xml
@@ -56,6 +56,11 @@ If the output is not quite correct, check for invisible
trailing spaces!
<release version="1.4" date="TBD" description="
New features, updates and bug fixes.
">
+ <action dev="aherbert" type="fix" due-to="Security scan, Alex Herbert"
issue="NUMBERS-209">
+ "Trigamma": Avoid an infinite loop on large negative arguments. All
negative
+ arguments are now computed using the reflection formula to map the
computation
+ to a positive argument.
+ </action>
</release>
<release version="1.3" date="2026-04-20" description="
