You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by er...@apache.org on 2016/05/14 16:23:11 UTC
[01/10] [math] Improved digamma function,
especially for negative values.
Repository: commons-math
Updated Branches:
refs/heads/develop f695c9ce3 -> 1325484c3
Improved digamma function, especially for negative values.
Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/53661585
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/53661585
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/53661585
Branch: refs/heads/develop
Commit: 5366158572bd88f18df7883ab823ea62b7701d77
Parents: cbae75b
Author: Eric Prescott-Gagnon <er...@jda.com>
Authored: Fri Jan 8 14:18:30 2016 -0500
Committer: Eric Prescott-Gagnon <er...@jda.com>
Committed: Thu May 5 13:48:13 2016 -0400
----------------------------------------------------------------------
.../org/apache/commons/math4/special/Gamma.java | 42 ++++++++++++--------
.../apache/commons/math4/special/GammaTest.java | 1 +
2 files changed, 26 insertions(+), 17 deletions(-)
----------------------------------------------------------------------
http://git-wip-us.apache.org/repos/asf/commons-math/blob/53661585/src/main/java/org/apache/commons/math4/special/Gamma.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math4/special/Gamma.java b/src/main/java/org/apache/commons/math4/special/Gamma.java
index 70aaec3..c59049c 100644
--- a/src/main/java/org/apache/commons/math4/special/Gamma.java
+++ b/src/main/java/org/apache/commons/math4/special/Gamma.java
@@ -427,18 +427,16 @@ public class Gamma {
* <p>Computes the digamma function of x.</p>
*
* <p>This is an independently written implementation of the algorithm described in
- * Jose Bernardo, Algorithm AS 103: Psi (Digamma) Function, Applied Statistics, 1976.</p>
+ * Jose Bernardo, Algorithm AS 103: Psi (Digamma) Function, Applied Statistics, 1976.
+ * A reflection formula (https://en.wikipedia.org/wiki/Digamma_function#Reflection_formula)
+ * is incorporated to improve performance on negative values.</p>
*
* <p>Some of the constants have been changed to increase accuracy at the moderate expense
- * of run-time. The result should be accurate to within 10^-8 absolute tolerance for
- * x >= 10^-5 and within 10^-8 relative tolerance for x > 0.</p>
- *
- * <p>Performance for large negative values of x will be quite expensive (proportional to
- * |x|). Accuracy for negative values of x should be about 10^-8 absolute for results
- * less than 10^5 and 10^-8 relative for results larger than that.</p>
+ * of run-time. The result should be accurate to within 10^-8 relative tolerance for
+ * 0 < x < 10^-5 and within 10^-8 absolute tolerance otherwise.</p>
*
* @param x Argument.
- * @return digamma(x) to within 10-8 relative or absolute error whichever is smaller.
+ * @return digamma(x) to within 10-8 relative or absolute error whichever is larger.
* @see <a href="http://en.wikipedia.org/wiki/Digamma_function">Digamma</a>
* @see <a href="http://www.uv.es/~bernardo/1976AppStatist.pdf">Bernardo's original article </a>
* @since 2.0
@@ -448,22 +446,32 @@ public class Gamma {
return x;
}
+ double digamma = 0.0;
+ if (x < 0) {
+ // use reflection formula to fall back into positive values
+ digamma -= FastMath.PI / FastMath.tan(FastMath.PI * x);
+ x = 1 - x;
+ }
+
if (x > 0 && x <= S_LIMIT) {
// use method 5 from Bernardo AS103
// accurate to O(x)
- return -GAMMA - 1 / x;
+ return digamma -GAMMA - 1 / x;
}
- if (x >= C_LIMIT) {
- // use method 4 (accurate to O(1/x^8)
- double inv = 1 / (x * x);
- // 1 1 1 1
- // log(x) - --- - ------ + ------- - -------
- // 2 x 12 x^2 120 x^4 252 x^6
- return FastMath.log(x) - 0.5 / x - inv * ((1.0 / 12) + inv * (1.0 / 120 - inv / 252));
+ while (x < C_LIMIT) {
+ digamma -= 1.0 / x;
+ x += 1.0;
}
- return digamma(x + 1) - 1 / x;
+ // use method 4 (accurate to O(1/x^8)
+ double inv = 1 / (x * x);
+ // 1 1 1 1
+ // log(x) - --- - ------ + ------- - -------
+ // 2 x 12 x^2 120 x^4 252 x^6
+ digamma += FastMath.log(x) - 0.5 / x - inv * ((1.0 / 12) + inv * (1.0 / 120 - inv / 252));
+
+ return digamma;
}
/**
http://git-wip-us.apache.org/repos/asf/commons-math/blob/53661585/src/test/java/org/apache/commons/math4/special/GammaTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math4/special/GammaTest.java b/src/test/java/org/apache/commons/math4/special/GammaTest.java
index 51982a3..c320526 100644
--- a/src/test/java/org/apache/commons/math4/special/GammaTest.java
+++ b/src/test/java/org/apache/commons/math4/special/GammaTest.java
@@ -107,6 +107,7 @@ public class GammaTest {
Assert.assertEquals(-100.56088545786867450, Gamma.digamma(0.01), eps);
Assert.assertEquals(-4.0390398965921882955, Gamma.digamma(-0.8), eps);
Assert.assertEquals(4.2003210041401844726, Gamma.digamma(-6.3), eps);
+ Assert.assertEquals(-3.110625123035E-5, Gamma.digamma(1.4616), eps);
}
@Test
[09/10] [math] Merge branch 'feature-MATH-1318' into develop
Posted by er...@apache.org.
Merge branch 'feature-MATH-1318' into develop
Completes issue MATH-1318 (see JIRA).
Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/64b47d86
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/64b47d86
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/64b47d86
Branch: refs/heads/develop
Commit: 64b47d86facb0dc2355a6295f143410938740ad3
Parents: fb1c112 db50bb3
Author: Gilles <er...@apache.org>
Authored: Sat May 14 16:26:21 2016 +0200
Committer: Gilles <er...@apache.org>
Committed: Sat May 14 16:26:21 2016 +0200
----------------------------------------------------------------------
src/changes/changes.xml | 4 ++
.../org/apache/commons/math4/special/Gamma.java | 44 ++++++++++++--------
.../apache/commons/math4/special/GammaTest.java | 1 +
3 files changed, 31 insertions(+), 18 deletions(-)
----------------------------------------------------------------------
http://git-wip-us.apache.org/repos/asf/commons-math/blob/64b47d86/src/changes/changes.xml
----------------------------------------------------------------------
diff --cc src/changes/changes.xml
index 0bc6938,b33777d..897ca51
--- a/src/changes/changes.xml
+++ b/src/changes/changes.xml
@@@ -54,9 -54,10 +54,13 @@@ If the output is not quite correct, che
</release>
<release version="4.0" date="XXXX-XX-XX" description="">
+ <action dev="erans" type="add" issue="MATH-1015" due-to="Thomas Neidhart">
+ Gauss-Laguerre quadrature.
+ </action>
+ <action dev="erans" type="update" issue="MATH-1318" due-to="Eric Prescott-Gagnon">
+ "o.a.c.m.special.Gamma.digamma": Improved performance (through the use of
+ the reflection formula for negative arguments).
+ </action>
<action dev="erans" type="add" issue="MATH-1350" due-to="Rob Tompkins">
Improved code coverage (unit tests).
</action>
[03/10] [math] MATH-1015
Posted by er...@apache.org.
MATH-1015
Gauss-Laguerre quadrature scheme.
Thanks to Thomas Neidhart.
Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/93d35c8f
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/93d35c8f
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/93d35c8f
Branch: refs/heads/develop
Commit: 93d35c8f291b8953f8460c0eb02a51cc6d6cac22
Parents: cbae75b
Author: Gilles <er...@apache.org>
Authored: Sat May 7 01:25:42 2016 +0200
Committer: Gilles <er...@apache.org>
Committed: Sat May 7 01:25:42 2016 +0200
----------------------------------------------------------------------
.../gauss/GaussIntegratorFactory.java | 21 +++++
.../integration/gauss/LaguerreRuleFactory.java | 84 ++++++++++++++++++++
.../integration/gauss/LaguerreTest.java | 50 ++++++++++++
3 files changed, 155 insertions(+)
----------------------------------------------------------------------
http://git-wip-us.apache.org/repos/asf/commons-math/blob/93d35c8f/src/main/java/org/apache/commons/math4/analysis/integration/gauss/GaussIntegratorFactory.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math4/analysis/integration/gauss/GaussIntegratorFactory.java b/src/main/java/org/apache/commons/math4/analysis/integration/gauss/GaussIntegratorFactory.java
index d15d6ab..4dacc67 100644
--- a/src/main/java/org/apache/commons/math4/analysis/integration/gauss/GaussIntegratorFactory.java
+++ b/src/main/java/org/apache/commons/math4/analysis/integration/gauss/GaussIntegratorFactory.java
@@ -35,6 +35,27 @@ public class GaussIntegratorFactory {
private final BaseRuleFactory<BigDecimal> legendreHighPrecision = new LegendreHighPrecisionRuleFactory();
/** Generator of Gauss-Hermite integrators. */
private final BaseRuleFactory<Double> hermite = new HermiteRuleFactory();
+ /** Generator of Gauss-Laguerre integrators. */
+ private final BaseRuleFactory<Double> laguerre = new LaguerreRuleFactory();
+
+ /**
+ * Creates a Gauss-Laguerre integrator of the given order.
+ * The call to the
+ * {@link GaussIntegrator#integrate(org.apache.commons.math4.analysis.UnivariateFunction)
+ * integrate} method will perform an integration on the interval
+ * \([0, +\infty)\): the computed value is the improper integral of
+ * \(e^{-x} f(x)\)
+ * where \(f(x)\) is the function passed to the
+ * {@link SymmetricGaussIntegrator#integrate(org.apache.commons.math4.analysis.UnivariateFunction)
+ * integrate} method.
+ *
+ * @param numberOfPoints Order of the integration rule.
+ * @return a Gauss-Legendre integrator.
+ * @since 4.0
+ */
+ public GaussIntegrator laguerre(int numberOfPoints) {
+ return new GaussIntegrator(getRule(laguerre, numberOfPoints));
+ }
/**
* Creates a Gauss-Legendre integrator of the given order.
http://git-wip-us.apache.org/repos/asf/commons-math/blob/93d35c8f/src/main/java/org/apache/commons/math4/analysis/integration/gauss/LaguerreRuleFactory.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math4/analysis/integration/gauss/LaguerreRuleFactory.java b/src/main/java/org/apache/commons/math4/analysis/integration/gauss/LaguerreRuleFactory.java
new file mode 100644
index 0000000..220c460
--- /dev/null
+++ b/src/main/java/org/apache/commons/math4/analysis/integration/gauss/LaguerreRuleFactory.java
@@ -0,0 +1,84 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math4.analysis.integration.gauss;
+
+import java.util.Arrays;
+
+import org.apache.commons.math4.analysis.polynomials.PolynomialFunction;
+import org.apache.commons.math4.analysis.polynomials.PolynomialsUtils;
+import org.apache.commons.math4.exception.DimensionMismatchException;
+import org.apache.commons.math4.linear.EigenDecomposition;
+import org.apache.commons.math4.linear.MatrixUtils;
+import org.apache.commons.math4.linear.RealMatrix;
+import org.apache.commons.math4.util.Pair;
+
+/**
+ * Factory that creates Gauss-type quadrature rule using Laguerre polynomials.
+ *
+ * @see <a href="http://en.wikipedia.org/wiki/Gauss%E2%80%93Laguerre_quadrature">Gauss-Laguerre quadrature (Wikipedia)</a>
+ * @since 4.0
+ */
+public class LaguerreRuleFactory extends BaseRuleFactory<Double> {
+ /** {@inheritDoc} */
+ @Override
+ protected Pair<Double[], Double[]> computeRule(int numberOfPoints)
+ throws DimensionMismatchException {
+
+ final RealMatrix companionMatrix = companionMatrix(numberOfPoints);
+ final EigenDecomposition eigen = new EigenDecomposition(companionMatrix);
+ final double[] roots = eigen.getRealEigenvalues();
+ Arrays.sort(roots);
+
+ final Double[] points = new Double[numberOfPoints];
+ final Double[] weights = new Double[numberOfPoints];
+
+ final int n1 = numberOfPoints + 1;
+ final long n1Squared = n1 * (long) n1;
+ final PolynomialFunction laguerreN1 = PolynomialsUtils.createLaguerrePolynomial(n1);
+ for (int i = 0; i < numberOfPoints; i++) {
+ final double xi = roots[i];
+ points[i] = xi;
+
+ final double val = laguerreN1.value(xi);
+ weights[i] = xi / n1Squared / (val * val);
+ }
+
+ return new Pair<Double[], Double[]>(points, weights);
+ }
+
+ /**
+ * @param degree Matrix dimension.
+ * @return a square matrix.
+ */
+ private RealMatrix companionMatrix(final int degree) {
+ final RealMatrix c = MatrixUtils.createRealMatrix(degree, degree);
+
+ for (int i = 0; i < degree; i++) {
+ c.setEntry(i, i, 2 * i + 1);
+ if (i + 1 < degree) {
+ // subdiagonal
+ c.setEntry(i+1, i, -(i + 1));
+ }
+ if (i - 1 >= 0) {
+ // superdiagonal
+ c.setEntry(i-1, i, -i);
+ }
+ }
+
+ return c;
+ }
+}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/93d35c8f/src/test/java/org/apache/commons/math4/analysis/integration/gauss/LaguerreTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math4/analysis/integration/gauss/LaguerreTest.java b/src/test/java/org/apache/commons/math4/analysis/integration/gauss/LaguerreTest.java
new file mode 100644
index 0000000..87b502b
--- /dev/null
+++ b/src/test/java/org/apache/commons/math4/analysis/integration/gauss/LaguerreTest.java
@@ -0,0 +1,50 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math4.analysis.integration.gauss;
+
+import org.apache.commons.math4.analysis.UnivariateFunction;
+import org.apache.commons.math4.special.Gamma;
+import org.apache.commons.math4.util.FastMath;
+import org.junit.Assert;
+import org.junit.Test;
+
+/**
+ * Test of the {@link LaguerreRuleFactory}.
+ */
+public class LaguerreTest {
+ private static final GaussIntegratorFactory factory = new GaussIntegratorFactory();
+
+ @Test
+ public void testGamma() {
+ final double tol = 1e-13;
+
+ for (int i = 2; i < 10; i += 1) {
+ final double t = i;
+
+ final UnivariateFunction f = new UnivariateFunction() {
+ @Override
+ public double value(double x) {
+ return FastMath.pow(x, t - 1);
+ }
+ };
+
+ final GaussIntegrator integrator = factory.laguerre(7);
+ final double s = integrator.integrate(f);
+ Assert.assertEquals(1d, Gamma.gamma(t) / s, tol);
+ }
+ }
+}
[06/10] [math] Update "changes.xml".
Posted by er...@apache.org.
Update "changes.xml".
Attribution.
Trailing spaces removed.
Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/aa892624
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/aa892624
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/aa892624
Branch: refs/heads/develop
Commit: aa8926241048b4d6cbfc981feebb5d19d5799e3d
Parents: b77184a
Author: Gilles <er...@apache.org>
Authored: Tue May 10 23:52:51 2016 +0200
Committer: Gilles <er...@apache.org>
Committed: Tue May 10 23:52:51 2016 +0200
----------------------------------------------------------------------
src/changes/changes.xml | 15 +++++++++------
1 file changed, 9 insertions(+), 6 deletions(-)
----------------------------------------------------------------------
http://git-wip-us.apache.org/repos/asf/commons-math/blob/aa892624/src/changes/changes.xml
----------------------------------------------------------------------
diff --git a/src/changes/changes.xml b/src/changes/changes.xml
index f8b5658..0bc6938 100644
--- a/src/changes/changes.xml
+++ b/src/changes/changes.xml
@@ -54,6 +54,9 @@ If the output is not quite correct, check for invisible trailing spaces!
</release>
<release version="4.0" date="XXXX-XX-XX" description="">
+ <action dev="erans" type="add" issue="MATH-1015" due-to="Thomas Neidhart">
+ Gauss-Laguerre quadrature.
+ </action>
<action dev="erans" type="add" issue="MATH-1350" due-to="Rob Tompkins">
Improved code coverage (unit tests).
</action>
@@ -141,35 +144,35 @@ If the output is not quite correct, check for invisible trailing spaces!
"MathRuntimeException" is now the base class for all commons-math
exceptions (except for "NullArgumentException" which extends
"NullPointerException").
- </action>
+ </action>
<action dev="tn" type="remove" issue="MATH-1205">
Removed methods "test(...)" from "AbstractUnivariateStatistic".
The already existing methods "MathArrays#verifyValues(...)" shall
be used instead.
- </action>
+ </action>
<action dev="tn" type="update" issue="MATH-1205">
The abstract class "AbstractStorelessUnivariateStatistic" does not
extend anymore from "AbstractUnivariateStatistic".
- </action>
+ </action>
<action dev="tn" type="update" issue="MATH-1205">
Default implementation of
"AbstractStorelessUnivariateStatistic#equals(Object)"
will only return true if both instances have the same type. Previously
different statistics were considered to be equal if their current state
happened to be equal.
- </action>
+ </action>
<action dev="tn" type="update" issue="MATH-1205">
Default implementations of "AbstractStorelessUnivariateStatistic#evaluate(...)"
do not alter the internal state anymore. Instead a temporary copy of
the statistic is created for evaluation purposes.
- </action>
+ </action>
<action dev="tn" type="fix" issue="MATH-1205">
Methods "evaluate(...)" of class "Variance" changed the internal state
although it was stated differently in the javadoc.
</action>
<action dev="luc" type="fix" issue="MATH-1191">
Fixed ignored method parameters in QRDecomposition protected methods.
- </action>
+ </action>
<action dev="luc" type="fix" issue="MATH-1212">
Changed javadoc as the RandomDataGenerator class does not implement
an interface anymore (the previous interface has been deprecated in
[05/10] [math] MATH-1015
Posted by er...@apache.org.
MATH-1015
Gauss-Laguerre quadrature scheme.
Thanks to Thomas Neidhart.
Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/b77184a2
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/b77184a2
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/b77184a2
Branch: refs/heads/develop
Commit: b77184a2189c8f91100d148354d1fdf38e83d13f
Parents: f695c9c
Author: Gilles <er...@apache.org>
Authored: Sat May 7 01:25:42 2016 +0200
Committer: Gilles <er...@apache.org>
Committed: Tue May 10 23:47:36 2016 +0200
----------------------------------------------------------------------
.../gauss/GaussIntegratorFactory.java | 21 +++++
.../integration/gauss/LaguerreRuleFactory.java | 84 ++++++++++++++++++++
.../integration/gauss/LaguerreTest.java | 50 ++++++++++++
3 files changed, 155 insertions(+)
----------------------------------------------------------------------
http://git-wip-us.apache.org/repos/asf/commons-math/blob/b77184a2/src/main/java/org/apache/commons/math4/analysis/integration/gauss/GaussIntegratorFactory.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math4/analysis/integration/gauss/GaussIntegratorFactory.java b/src/main/java/org/apache/commons/math4/analysis/integration/gauss/GaussIntegratorFactory.java
index d15d6ab..4dacc67 100644
--- a/src/main/java/org/apache/commons/math4/analysis/integration/gauss/GaussIntegratorFactory.java
+++ b/src/main/java/org/apache/commons/math4/analysis/integration/gauss/GaussIntegratorFactory.java
@@ -35,6 +35,27 @@ public class GaussIntegratorFactory {
private final BaseRuleFactory<BigDecimal> legendreHighPrecision = new LegendreHighPrecisionRuleFactory();
/** Generator of Gauss-Hermite integrators. */
private final BaseRuleFactory<Double> hermite = new HermiteRuleFactory();
+ /** Generator of Gauss-Laguerre integrators. */
+ private final BaseRuleFactory<Double> laguerre = new LaguerreRuleFactory();
+
+ /**
+ * Creates a Gauss-Laguerre integrator of the given order.
+ * The call to the
+ * {@link GaussIntegrator#integrate(org.apache.commons.math4.analysis.UnivariateFunction)
+ * integrate} method will perform an integration on the interval
+ * \([0, +\infty)\): the computed value is the improper integral of
+ * \(e^{-x} f(x)\)
+ * where \(f(x)\) is the function passed to the
+ * {@link SymmetricGaussIntegrator#integrate(org.apache.commons.math4.analysis.UnivariateFunction)
+ * integrate} method.
+ *
+ * @param numberOfPoints Order of the integration rule.
+ * @return a Gauss-Legendre integrator.
+ * @since 4.0
+ */
+ public GaussIntegrator laguerre(int numberOfPoints) {
+ return new GaussIntegrator(getRule(laguerre, numberOfPoints));
+ }
/**
* Creates a Gauss-Legendre integrator of the given order.
http://git-wip-us.apache.org/repos/asf/commons-math/blob/b77184a2/src/main/java/org/apache/commons/math4/analysis/integration/gauss/LaguerreRuleFactory.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math4/analysis/integration/gauss/LaguerreRuleFactory.java b/src/main/java/org/apache/commons/math4/analysis/integration/gauss/LaguerreRuleFactory.java
new file mode 100644
index 0000000..220c460
--- /dev/null
+++ b/src/main/java/org/apache/commons/math4/analysis/integration/gauss/LaguerreRuleFactory.java
@@ -0,0 +1,84 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math4.analysis.integration.gauss;
+
+import java.util.Arrays;
+
+import org.apache.commons.math4.analysis.polynomials.PolynomialFunction;
+import org.apache.commons.math4.analysis.polynomials.PolynomialsUtils;
+import org.apache.commons.math4.exception.DimensionMismatchException;
+import org.apache.commons.math4.linear.EigenDecomposition;
+import org.apache.commons.math4.linear.MatrixUtils;
+import org.apache.commons.math4.linear.RealMatrix;
+import org.apache.commons.math4.util.Pair;
+
+/**
+ * Factory that creates Gauss-type quadrature rule using Laguerre polynomials.
+ *
+ * @see <a href="http://en.wikipedia.org/wiki/Gauss%E2%80%93Laguerre_quadrature">Gauss-Laguerre quadrature (Wikipedia)</a>
+ * @since 4.0
+ */
+public class LaguerreRuleFactory extends BaseRuleFactory<Double> {
+ /** {@inheritDoc} */
+ @Override
+ protected Pair<Double[], Double[]> computeRule(int numberOfPoints)
+ throws DimensionMismatchException {
+
+ final RealMatrix companionMatrix = companionMatrix(numberOfPoints);
+ final EigenDecomposition eigen = new EigenDecomposition(companionMatrix);
+ final double[] roots = eigen.getRealEigenvalues();
+ Arrays.sort(roots);
+
+ final Double[] points = new Double[numberOfPoints];
+ final Double[] weights = new Double[numberOfPoints];
+
+ final int n1 = numberOfPoints + 1;
+ final long n1Squared = n1 * (long) n1;
+ final PolynomialFunction laguerreN1 = PolynomialsUtils.createLaguerrePolynomial(n1);
+ for (int i = 0; i < numberOfPoints; i++) {
+ final double xi = roots[i];
+ points[i] = xi;
+
+ final double val = laguerreN1.value(xi);
+ weights[i] = xi / n1Squared / (val * val);
+ }
+
+ return new Pair<Double[], Double[]>(points, weights);
+ }
+
+ /**
+ * @param degree Matrix dimension.
+ * @return a square matrix.
+ */
+ private RealMatrix companionMatrix(final int degree) {
+ final RealMatrix c = MatrixUtils.createRealMatrix(degree, degree);
+
+ for (int i = 0; i < degree; i++) {
+ c.setEntry(i, i, 2 * i + 1);
+ if (i + 1 < degree) {
+ // subdiagonal
+ c.setEntry(i+1, i, -(i + 1));
+ }
+ if (i - 1 >= 0) {
+ // superdiagonal
+ c.setEntry(i-1, i, -i);
+ }
+ }
+
+ return c;
+ }
+}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/b77184a2/src/test/java/org/apache/commons/math4/analysis/integration/gauss/LaguerreTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math4/analysis/integration/gauss/LaguerreTest.java b/src/test/java/org/apache/commons/math4/analysis/integration/gauss/LaguerreTest.java
new file mode 100644
index 0000000..87b502b
--- /dev/null
+++ b/src/test/java/org/apache/commons/math4/analysis/integration/gauss/LaguerreTest.java
@@ -0,0 +1,50 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math4.analysis.integration.gauss;
+
+import org.apache.commons.math4.analysis.UnivariateFunction;
+import org.apache.commons.math4.special.Gamma;
+import org.apache.commons.math4.util.FastMath;
+import org.junit.Assert;
+import org.junit.Test;
+
+/**
+ * Test of the {@link LaguerreRuleFactory}.
+ */
+public class LaguerreTest {
+ private static final GaussIntegratorFactory factory = new GaussIntegratorFactory();
+
+ @Test
+ public void testGamma() {
+ final double tol = 1e-13;
+
+ for (int i = 2; i < 10; i += 1) {
+ final double t = i;
+
+ final UnivariateFunction f = new UnivariateFunction() {
+ @Override
+ public double value(double x) {
+ return FastMath.pow(x, t - 1);
+ }
+ };
+
+ final GaussIntegrator integrator = factory.laguerre(7);
+ final double s = integrator.integrate(f);
+ Assert.assertEquals(1d, Gamma.gamma(t) / s, tol);
+ }
+ }
+}
[08/10] [math] Merge branch 'feature-MATH-1015' into develop
Posted by er...@apache.org.
Merge branch 'feature-MATH-1015' into develop
Completes issue MATH-1015 (see JIRA).
Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/fb1c1127
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/fb1c1127
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/fb1c1127
Branch: refs/heads/develop
Commit: fb1c11278810919fb5cbbadd5e7c3167cd4fc055
Parents: f695c9c bbd702a
Author: Gilles <er...@apache.org>
Authored: Sat May 14 15:17:40 2016 +0200
Committer: Gilles <er...@apache.org>
Committed: Sat May 14 15:17:40 2016 +0200
----------------------------------------------------------------------
src/changes/changes.xml | 15 ++--
.../gauss/GaussIntegratorFactory.java | 21 +++++
.../integration/gauss/LaguerreRuleFactory.java | 84 ++++++++++++++++++++
.../integration/gauss/LaguerreTest.java | 50 ++++++++++++
4 files changed, 164 insertions(+), 6 deletions(-)
----------------------------------------------------------------------
[02/10] [math] Update "changes.xml".
Posted by er...@apache.org.
Update "changes.xml".
Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/9e877986
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/9e877986
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/9e877986
Branch: refs/heads/develop
Commit: 9e877986cd976df2ec84118c523174c0f35ef89d
Parents: 5366158
Author: Gilles <er...@apache.org>
Authored: Fri May 6 12:36:38 2016 +0200
Committer: Gilles <er...@apache.org>
Committed: Fri May 6 12:36:38 2016 +0200
----------------------------------------------------------------------
src/changes/changes.xml | 4 ++++
1 file changed, 4 insertions(+)
----------------------------------------------------------------------
http://git-wip-us.apache.org/repos/asf/commons-math/blob/9e877986/src/changes/changes.xml
----------------------------------------------------------------------
diff --git a/src/changes/changes.xml b/src/changes/changes.xml
index f8b5658..b33777d 100644
--- a/src/changes/changes.xml
+++ b/src/changes/changes.xml
@@ -54,6 +54,10 @@ If the output is not quite correct, check for invisible trailing spaces!
</release>
<release version="4.0" date="XXXX-XX-XX" description="">
+ <action dev="erans" type="update" issue="MATH-1318" due-to="Eric Prescott-Gagnon">
+ "o.a.c.m.special.Gamma.digamma": Improved performance (through the use of
+ the reflection formula for negative arguments).
+ </action>
<action dev="erans" type="add" issue="MATH-1350" due-to="Rob Tompkins">
Improved code coverage (unit tests).
</action>
[10/10] [math] Update "changes.xml".
Posted by er...@apache.org.
Update "changes.xml".
Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/1325484c
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/1325484c
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/1325484c
Branch: refs/heads/develop
Commit: 1325484c33b933f7da1cbc1ac6e289ee65a61511
Parents: 64b47d8
Author: Gilles <er...@apache.org>
Authored: Sat May 14 16:52:02 2016 +0200
Committer: Gilles <er...@apache.org>
Committed: Sat May 14 16:52:02 2016 +0200
----------------------------------------------------------------------
src/changes/changes.xml | 5 +++++
1 file changed, 5 insertions(+)
----------------------------------------------------------------------
http://git-wip-us.apache.org/repos/asf/commons-math/blob/1325484c/src/changes/changes.xml
----------------------------------------------------------------------
diff --git a/src/changes/changes.xml b/src/changes/changes.xml
index 897ca51..3ced6eb 100644
--- a/src/changes/changes.xml
+++ b/src/changes/changes.xml
@@ -64,6 +64,11 @@ If the output is not quite correct, check for invisible trailing spaces!
<action dev="erans" type="add" issue="MATH-1350" due-to="Rob Tompkins">
Improved code coverage (unit tests).
</action>
+ <action dev="erans" type="add" issue="MATH-1336">
+ New 64-bits RNG implementations. On 64-bits systems, they are ~35% faster than
+ "java.util.Random" for generating "long" or "double" values. They also guarantee
+ much better randomness than "Random" (cf. MATH-1327).
+ </action>
<action dev="erans" type="add" issue="MATH-1335">
Refactoring of uniform random number generator functionality: new API
implemented in package "o.a.c.m.rng".
[07/10] [math] Merge branch 'feature-MATH-1015' of
https://git-wip-us.apache.org/repos/asf/commons-math into feature-MATH-1015
Posted by er...@apache.org.
Merge branch 'feature-MATH-1015' of https://git-wip-us.apache.org/repos/asf/commons-math into feature-MATH-1015
Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/bbd702ae
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/bbd702ae
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/bbd702ae
Branch: refs/heads/develop
Commit: bbd702ae2f6c52b58d2e2ba99d7fc548545f46e6
Parents: aa89262 93d35c8
Author: Gilles <er...@apache.org>
Authored: Tue May 10 23:55:01 2016 +0200
Committer: Gilles <er...@apache.org>
Committed: Tue May 10 23:55:01 2016 +0200
----------------------------------------------------------------------
----------------------------------------------------------------------
[04/10] [math] Typo.
Posted by er...@apache.org.
Typo.
Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/db50bb3c
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/db50bb3c
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/db50bb3c
Branch: refs/heads/develop
Commit: db50bb3c3c3311bf286dcebf1ffe6edef25f2aee
Parents: 9e87798
Author: Gilles <er...@apache.org>
Authored: Tue May 10 23:41:59 2016 +0200
Committer: Gilles <er...@apache.org>
Committed: Tue May 10 23:41:59 2016 +0200
----------------------------------------------------------------------
src/main/java/org/apache/commons/math4/special/Gamma.java | 4 ++--
1 file changed, 2 insertions(+), 2 deletions(-)
----------------------------------------------------------------------
http://git-wip-us.apache.org/repos/asf/commons-math/blob/db50bb3c/src/main/java/org/apache/commons/math4/special/Gamma.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math4/special/Gamma.java b/src/main/java/org/apache/commons/math4/special/Gamma.java
index c59049c..e798277 100644
--- a/src/main/java/org/apache/commons/math4/special/Gamma.java
+++ b/src/main/java/org/apache/commons/math4/special/Gamma.java
@@ -436,7 +436,7 @@ public class Gamma {
* 0 < x < 10^-5 and within 10^-8 absolute tolerance otherwise.</p>
*
* @param x Argument.
- * @return digamma(x) to within 10-8 relative or absolute error whichever is larger.
+ * @return digamma(x) to within 10^-8 relative or absolute error whichever is larger.
* @see <a href="http://en.wikipedia.org/wiki/Digamma_function">Digamma</a>
* @see <a href="http://www.uv.es/~bernardo/1976AppStatist.pdf">Bernardo's original article </a>
* @since 2.0
@@ -480,7 +480,7 @@ public class Gamma {
* of digamma.
*
* @param x Argument.
- * @return trigamma(x) to within 10-8 relative or absolute error whichever is smaller
+ * @return trigamma(x) to within 10^-8 relative or absolute error whichever is smaller
* @see <a href="http://en.wikipedia.org/wiki/Trigamma_function">Trigamma</a>
* @see Gamma#digamma(double)
* @since 2.0