You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@mahout.apache.org by ra...@apache.org on 2018/09/08 23:35:10 UTC
[06/15] mahout git commit: NO-JIRA Trevors updates
http://git-wip-us.apache.org/repos/asf/mahout/blob/545648f6/core/src/main/java/org/apache/mahout/math/jet/math/Polynomial.java
----------------------------------------------------------------------
diff --git a/core/src/main/java/org/apache/mahout/math/jet/math/Polynomial.java b/core/src/main/java/org/apache/mahout/math/jet/math/Polynomial.java
new file mode 100644
index 0000000..723e7d0
--- /dev/null
+++ b/core/src/main/java/org/apache/mahout/math/jet/math/Polynomial.java
@@ -0,0 +1,98 @@
+/*
+ * 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.
+ */
+
+/*
+Copyright 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
+is hereby granted without fee, provided that the above copyright notice appear in all copies and
+that both that copyright notice and this permission notice appear in supporting documentation.
+CERN makes no representations about the suitability of this software for any purpose.
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.math.jet.math;
+
+/**
+ * Polynomial functions.
+ */
+public final class Polynomial {
+
+ private Polynomial() {
+ }
+
+ /**
+ * Evaluates the given polynomial of degree <tt>N</tt> at <tt>x</tt>, assuming coefficient of N is 1.0. Otherwise same
+ * as <tt>polevl()</tt>.
+ * <pre>
+ * 2 N
+ * y = C + C x + C x +...+ C x
+ * 0 1 2 N
+ *
+ * where C = 1 and hence is omitted from the array.
+ * N
+ *
+ * Coefficients are stored in reverse order:
+ *
+ * coef[0] = C , ..., coef[N-1] = C .
+ * N-1 0
+ *
+ * Calling arguments are otherwise the same as polevl().
+ * </pre>
+ * In the interest of speed, there are no checks for out of bounds arithmetic.
+ *
+ * @param x argument to the polynomial.
+ * @param coef the coefficients of the polynomial.
+ * @param N the degree of the polynomial.
+ */
+ public static double p1evl(double x, double[] coef, int N) {
+
+ double ans = x + coef[0];
+
+ for (int i = 1; i < N; i++) {
+ ans = ans * x + coef[i];
+ }
+
+ return ans;
+ }
+
+ /**
+ * Evaluates the given polynomial of degree <tt>N</tt> at <tt>x</tt>.
+ * <pre>
+ * 2 N
+ * y = C + C x + C x +...+ C x
+ * 0 1 2 N
+ *
+ * Coefficients are stored in reverse order:
+ *
+ * coef[0] = C , ..., coef[N] = C .
+ * N 0
+ * </pre>
+ * In the interest of speed, there are no checks for out of bounds arithmetic.
+ *
+ * @param x argument to the polynomial.
+ * @param coef the coefficients of the polynomial.
+ * @param N the degree of the polynomial.
+ */
+ public static double polevl(double x, double[] coef, int N) {
+ double ans = coef[0];
+
+ for (int i = 1; i <= N; i++) {
+ ans = ans * x + coef[i];
+ }
+
+ return ans;
+ }
+}
http://git-wip-us.apache.org/repos/asf/mahout/blob/545648f6/core/src/main/java/org/apache/mahout/math/jet/math/package-info.java
----------------------------------------------------------------------
diff --git a/core/src/main/java/org/apache/mahout/math/jet/math/package-info.java b/core/src/main/java/org/apache/mahout/math/jet/math/package-info.java
new file mode 100644
index 0000000..3cda850
--- /dev/null
+++ b/core/src/main/java/org/apache/mahout/math/jet/math/package-info.java
@@ -0,0 +1,5 @@
+/**
+ * Tools for basic and advanced mathematics: Arithmetics and Algebra, Polynomials and Chebyshev series, Bessel and Airy
+ * functions, Function Objects for generic function evaluation, etc.
+ */
+package org.apache.mahout.math.jet.math;
http://git-wip-us.apache.org/repos/asf/mahout/blob/545648f6/core/src/main/java/org/apache/mahout/math/jet/random/AbstractContinousDistribution.java
----------------------------------------------------------------------
diff --git a/core/src/main/java/org/apache/mahout/math/jet/random/AbstractContinousDistribution.java b/core/src/main/java/org/apache/mahout/math/jet/random/AbstractContinousDistribution.java
new file mode 100644
index 0000000..8ca03d0
--- /dev/null
+++ b/core/src/main/java/org/apache/mahout/math/jet/random/AbstractContinousDistribution.java
@@ -0,0 +1,51 @@
+/**
+ * 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.
+ */
+/*
+Copyright 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
+is hereby granted without fee, provided that the above copyright notice appear in all copies and
+that both that copyright notice and this permission notice appear in supporting documentation.
+CERN makes no representations about the suitability of this software for any purpose.
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.math.jet.random;
+
+/**
+ * Abstract base class for all continuous distributions. Continuous distributions have
+ * probability density and a cumulative distribution functions.
+ *
+ */
+public abstract class AbstractContinousDistribution extends AbstractDistribution {
+ public double cdf(double x) {
+ throw new UnsupportedOperationException("Can't compute pdf for " + this.getClass().getName());
+ }
+
+ public double pdf(double x) {
+ throw new UnsupportedOperationException("Can't compute pdf for " + this.getClass().getName());
+ }
+
+ /**
+ * @return A random number from the distribution; returns <tt>(int) Math.round(nextDouble())</tt>.
+ * Override this method if necessary.
+ */
+ @Override
+ public int nextInt() {
+ return (int) Math.round(nextDouble());
+ }
+}
http://git-wip-us.apache.org/repos/asf/mahout/blob/545648f6/core/src/main/java/org/apache/mahout/math/jet/random/AbstractDiscreteDistribution.java
----------------------------------------------------------------------
diff --git a/core/src/main/java/org/apache/mahout/math/jet/random/AbstractDiscreteDistribution.java b/core/src/main/java/org/apache/mahout/math/jet/random/AbstractDiscreteDistribution.java
new file mode 100644
index 0000000..d93d76c
--- /dev/null
+++ b/core/src/main/java/org/apache/mahout/math/jet/random/AbstractDiscreteDistribution.java
@@ -0,0 +1,27 @@
+/*
+Copyright 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
+is hereby granted without fee, provided that the above copyright notice appear in all copies and
+that both that copyright notice and this permission notice appear in supporting documentation.
+CERN makes no representations about the suitability of this software for any purpose.
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.math.jet.random;
+
+/**
+ * Abstract base class for all discrete distributions.
+ *
+ */
+public abstract class AbstractDiscreteDistribution extends AbstractDistribution {
+
+ /** Makes this class non instantiable, but still let's others inherit from it. */
+ protected AbstractDiscreteDistribution() {
+ }
+
+ /** Returns a random number from the distribution; returns <tt>(double) nextInt()</tt>. */
+ @Override
+ public double nextDouble() {
+ return nextInt();
+ }
+
+}
http://git-wip-us.apache.org/repos/asf/mahout/blob/545648f6/core/src/main/java/org/apache/mahout/math/jet/random/AbstractDistribution.java
----------------------------------------------------------------------
diff --git a/core/src/main/java/org/apache/mahout/math/jet/random/AbstractDistribution.java b/core/src/main/java/org/apache/mahout/math/jet/random/AbstractDistribution.java
new file mode 100644
index 0000000..8e9cb0e
--- /dev/null
+++ b/core/src/main/java/org/apache/mahout/math/jet/random/AbstractDistribution.java
@@ -0,0 +1,87 @@
+/*
+ * 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.
+ */
+
+/*
+Copyright � 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
+is hereby granted without fee, provided that the above copyright notice appear in all copies and
+that both that copyright notice and this permission notice appear in supporting documentation.
+CERN makes no representations about the suitability of this software for any purpose.
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.math.jet.random;
+
+import java.util.Random;
+
+import org.apache.mahout.math.function.DoubleFunction;
+import org.apache.mahout.math.function.IntFunction;
+
+public abstract class AbstractDistribution extends DoubleFunction implements IntFunction {
+
+ private Random randomGenerator;
+
+ /** Makes this class non instantiable, but still let's others inherit from it. */
+ protected AbstractDistribution() {
+ }
+
+ protected Random getRandomGenerator() {
+ return randomGenerator;
+ }
+
+ protected double randomDouble() {
+ return randomGenerator.nextDouble();
+ }
+
+ /**
+ * Equivalent to <tt>nextDouble()</tt>. This has the effect that distributions can now be used as function objects,
+ * returning a random number upon function evaluation.
+ */
+ @Override
+ public double apply(double dummy) {
+ return nextDouble();
+ }
+
+ /**
+ * Equivalent to <tt>nextInt()</tt>. This has the effect that distributions can now be used as function objects,
+ * returning a random number upon function evaluation.
+ */
+ @Override
+ public int apply(int dummy) {
+ return nextInt();
+ }
+
+ /**
+ * Returns a random number from the distribution.
+ * @return A new sample from this distribution.
+ */
+ public abstract double nextDouble();
+
+ /**
+ * @return
+ * A random number from the distribution; returns <tt>(int) Math.round(nextDouble())</tt>. Override this
+ * method if necessary.
+ */
+ public abstract int nextInt();
+
+ /**
+ * Sets the uniform random generator internally used.
+ * @param randomGenerator the new PRNG
+ */
+ public void setRandomGenerator(Random randomGenerator) {
+ this.randomGenerator = randomGenerator;
+ }
+}
http://git-wip-us.apache.org/repos/asf/mahout/blob/545648f6/core/src/main/java/org/apache/mahout/math/jet/random/Exponential.java
----------------------------------------------------------------------
diff --git a/core/src/main/java/org/apache/mahout/math/jet/random/Exponential.java b/core/src/main/java/org/apache/mahout/math/jet/random/Exponential.java
new file mode 100644
index 0000000..06472c2
--- /dev/null
+++ b/core/src/main/java/org/apache/mahout/math/jet/random/Exponential.java
@@ -0,0 +1,81 @@
+/*
+Copyright � 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
+is hereby granted without fee, provided that the above copyright notice appear in all copies and
+that both that copyright notice and this permission notice appear in supporting documentation.
+CERN makes no representations about the suitability of this software for any purpose.
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.math.jet.random;
+
+import java.util.Locale;
+import java.util.Random;
+
+public class Exponential extends AbstractContinousDistribution {
+ // rate parameter for the distribution. Mean is 1/lambda.
+ private double lambda;
+
+ /**
+ * Provides a negative exponential distribution given a rate parameter lambda and an underlying
+ * random number generator. The mean of this distribution will be equal to 1/lambda.
+ *
+ * @param lambda The rate parameter of the distribution.
+ * @param randomGenerator The PRNG that is used to generate values.
+ */
+ public Exponential(double lambda, Random randomGenerator) {
+ setRandomGenerator(randomGenerator);
+ this.lambda = lambda;
+ }
+
+ /**
+ * Returns the cumulative distribution function.
+ * @param x The point at which the cumulative distribution function is to be evaluated.
+ * @return Returns the integral from -infinity to x of the PDF, also known as the cumulative distribution
+ * function.
+ */
+ @Override
+ public double cdf(double x) {
+ if (x <= 0.0) {
+ return 0.0;
+ }
+ return 1.0 - Math.exp(-x * lambda);
+ }
+
+ /**
+ * Returns a random number from the distribution.
+ */
+ @Override
+ public double nextDouble() {
+ return -Math.log1p(-randomDouble()) / lambda;
+ }
+
+ /**
+ * Returns the value of the probability density function at a particular point.
+ * @param x The point at which the probability density function is to be evaluated.
+ * @return The value of the probability density function at the specified point.
+ */
+ @Override
+ public double pdf(double x) {
+ if (x < 0.0) {
+ return 0.0;
+ }
+ return lambda * Math.exp(-x * lambda);
+ }
+
+ /**
+ * Sets the rate parameter.
+ * @param lambda The new value of the rate parameter.
+ */
+ public void setState(double lambda) {
+ this.lambda = lambda;
+ }
+
+ /**
+ * Returns a String representation of the receiver.
+ */
+ @Override
+ public String toString() {
+ return String.format(Locale.ENGLISH, "%s(%.4f)", this.getClass().getName(), lambda);
+ }
+
+}
http://git-wip-us.apache.org/repos/asf/mahout/blob/545648f6/core/src/main/java/org/apache/mahout/math/jet/random/Gamma.java
----------------------------------------------------------------------
diff --git a/core/src/main/java/org/apache/mahout/math/jet/random/Gamma.java b/core/src/main/java/org/apache/mahout/math/jet/random/Gamma.java
new file mode 100644
index 0000000..914157b
--- /dev/null
+++ b/core/src/main/java/org/apache/mahout/math/jet/random/Gamma.java
@@ -0,0 +1,302 @@
+/*
+ * 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.
+ */
+
+/*
+Copyright � 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
+is hereby granted without fee, provided that the above copyright notice appear in all copies and
+that both that copyright notice and this permission notice appear in supporting documentation.
+CERN makes no representations about the suitability of this software for any purpose.
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.math.jet.random;
+
+import org.apache.mahout.math.jet.stat.Probability;
+
+import java.util.Random;
+
+public class Gamma extends AbstractContinousDistribution {
+ // shape
+ private final double alpha;
+
+ // rate
+ private final double rate;
+
+ /**
+ * Constructs a Gamma distribution with a given shape (alpha) and rate (beta).
+ *
+ * @param alpha The shape parameter.
+ * @param rate The rate parameter.
+ * @param randomGenerator The random number generator that generates bits for us.
+ * @throws IllegalArgumentException if <tt>alpha <= 0.0 || alpha <= 0.0</tt>.
+ */
+ public Gamma(double alpha, double rate, Random randomGenerator) {
+ this.alpha = alpha;
+ this.rate = rate;
+ setRandomGenerator(randomGenerator);
+ }
+
+ /**
+ * Returns the cumulative distribution function.
+ * @param x The end-point where the cumulation should end.
+ */
+ @Override
+ public double cdf(double x) {
+ return Probability.gamma(alpha, rate, x);
+ }
+
+ /** Returns a random number from the distribution. */
+ @Override
+ public double nextDouble() {
+ return nextDouble(alpha, rate);
+ }
+
+ /** Returns a random number from the distribution; bypasses the internal state.
+ * *
+ * Gamma Distribution - Acceptance Rejection combined with *
+ * Acceptance Complement *
+ * *
+ ******************************************************************
+ * *
+ * FUNCTION: - gds samples a random number from the standard *
+ * gamma distribution with parameter a > 0. *
+ * Acceptance Rejection gs for a < 1 , *
+ * Acceptance Complement gd for a >= 1 . *
+ * REFERENCES: - J.H. Ahrens, U. Dieter (1974): Computer methods *
+ * for sampling from gamma, beta, Poisson and *
+ * binomial distributions, Computing 12, 223-246. *
+ * - J.H. Ahrens, U. Dieter (1982): Generating gamma *
+ * variates by a modified rejection technique, *
+ * Communications of the ACM 25, 47-54. *
+ * SUBPROGRAMS: - drand(seed) ... (0,1)-Uniform generator with *
+ * unsigned long integer *seed *
+ * - NORMAL(seed) ... Normal generator N(0,1). *
+ * *
+ * @param alpha Shape parameter.
+ * @param rate Rate parameter (=1/scale).
+ * @return A gamma distributed sample.
+ */
+ public double nextDouble(double alpha, double rate) {
+ if (alpha <= 0.0) {
+ throw new IllegalArgumentException();
+ }
+ if (rate <= 0.0) {
+ throw new IllegalArgumentException();
+ }
+
+ double gds;
+ double b = 0.0;
+ if (alpha < 1.0) { // CASE A: Acceptance rejection algorithm gs
+ b = 1.0 + 0.36788794412 * alpha; // Step 1
+ while (true) {
+ double p = b * randomDouble();
+ if (p <= 1.0) { // Step 2. Case gds <= 1
+ gds = Math.exp(Math.log(p) / alpha);
+ if (Math.log(randomDouble()) <= -gds) {
+ return gds / rate;
+ }
+ } else { // Step 3. Case gds > 1
+ gds = -Math.log((b - p) / alpha);
+ if (Math.log(randomDouble()) <= (alpha - 1.0) * Math.log(gds)) {
+ return gds / rate;
+ }
+ }
+ }
+ } else { // CASE B: Acceptance complement algorithm gd (gaussian distribution, box muller transformation)
+ double ss = 0.0;
+ double s = 0.0;
+ double d = 0.0;
+ if (alpha != -1.0) { // Step 1. Preparations
+ ss = alpha - 0.5;
+ s = Math.sqrt(ss);
+ d = 5.656854249 - 12.0 * s;
+ }
+ // Step 2. Normal deviate
+ double v12;
+ double v1;
+ do {
+ v1 = 2.0 * randomDouble() - 1.0;
+ double v2 = 2.0 * randomDouble() - 1.0;
+ v12 = v1 * v1 + v2 * v2;
+ } while (v12 > 1.0);
+ double t = v1 * Math.sqrt(-2.0 * Math.log(v12) / v12);
+ double x = s + 0.5 * t;
+ gds = x * x;
+ if (t >= 0.0) {
+ return gds / rate;
+ } // Immediate acceptance
+
+ double u = randomDouble();
+ if (d * u <= t * t * t) {
+ return gds / rate;
+ } // Squeeze acceptance
+
+ double q0 = 0.0;
+ double si = 0.0;
+ double c = 0.0;
+ if (alpha != -1.0) { // Step 4. Set-up for hat case
+ double r = 1.0 / alpha;
+ double q9 = 0.0001710320;
+ double q8 = -0.0004701849;
+ double q7 = 0.0006053049;
+ double q6 = 0.0003340332;
+ double q5 = -0.0003349403;
+ double q4 = 0.0015746717;
+ double q3 = 0.0079849875;
+ double q2 = 0.0208333723;
+ double q1 = 0.0416666664;
+ q0 = ((((((((q9 * r + q8) * r + q7) * r + q6) * r + q5) * r + q4) * r + q3) * r + q2) * r + q1) * r;
+ if (alpha > 3.686) {
+ if (alpha > 13.022) {
+ b = 1.77;
+ si = 0.75;
+ c = 0.1515 / s;
+ } else {
+ b = 1.654 + 0.0076 * ss;
+ si = 1.68 / s + 0.275;
+ c = 0.062 / s + 0.024;
+ }
+ } else {
+ b = 0.463 + s - 0.178 * ss;
+ si = 1.235;
+ c = 0.195 / s - 0.079 + 0.016 * s;
+ }
+ }
+ double v;
+ double q;
+ double a9 = 0.104089866;
+ double a8 = -0.112750886;
+ double a7 = 0.110368310;
+ double a6 = -0.124385581;
+ double a5 = 0.142873973;
+ double a4 = -0.166677482;
+ double a3 = 0.199999867;
+ double a2 = -0.249999949;
+ double a1 = 0.333333333;
+ if (x > 0.0) { // Step 5. Calculation of q
+ v = t / (s + s); // Step 6.
+ if (Math.abs(v) > 0.25) {
+ q = q0 - s * t + 0.25 * t * t + (ss + ss) * Math.log1p(v);
+ } else {
+ q = q0 + 0.5 * t * t * ((((((((a9 * v + a8) * v + a7) * v + a6)
+ * v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v;
+ } // Step 7. Quotient acceptance
+ if (Math.log1p(-u) <= q) {
+ return gds / rate;
+ }
+ }
+
+ double e7 = 0.000247453;
+ double e6 = 0.001353826;
+ double e5 = 0.008345522;
+ double e4 = 0.041664508;
+ double e3 = 0.166666848;
+ double e2 = 0.499999994;
+ double e1 = 1.000000000;
+ while (true) { // Step 8. Double exponential deviate t
+ double sign_u;
+ double e;
+ do {
+ e = -Math.log(randomDouble());
+ u = randomDouble();
+ u = u + u - 1.0;
+ sign_u = u > 0 ? 1.0 : -1.0;
+ t = b + e * si * sign_u;
+ } while (t <= -0.71874483771719); // Step 9. Rejection of t
+ v = t / (s + s); // Step 10. New q(t)
+ if (Math.abs(v) > 0.25) {
+ q = q0 - s * t + 0.25 * t * t + (ss + ss) * Math.log1p(v);
+ } else {
+ q = q0 + 0.5 * t * t * ((((((((a9 * v + a8) * v + a7) * v + a6)
+ * v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v;
+ }
+ if (q <= 0.0) {
+ continue;
+ } // Step 11.
+ double w;
+ if (q > 0.5) {
+ w = Math.exp(q) - 1.0;
+ } else {
+ w = ((((((e7 * q + e6) * q + e5) * q + e4) * q + e3) * q + e2) * q + e1) * q;
+ } // Step 12. Hat acceptance
+ if (c * u * sign_u <= w * Math.exp(e - 0.5 * t * t)) {
+ x = s + 0.5 * t;
+ return x * x / rate;
+ }
+ }
+ }
+ }
+
+ /** Returns the probability distribution function.
+ * @param x Where to compute the density function.
+ *
+ * @return The value of the gamma density at x.
+ */
+ @Override
+ public double pdf(double x) {
+ if (x < 0) {
+ throw new IllegalArgumentException();
+ }
+ if (x == 0) {
+ if (alpha == 1.0) {
+ return rate;
+ } else if (alpha < 1) {
+ return Double.POSITIVE_INFINITY;
+ } else {
+ return 0;
+ }
+ }
+ if (alpha == 1.0) {
+ return rate * Math.exp(-x * rate);
+ }
+ return rate * Math.exp((alpha - 1.0) * Math.log(x * rate) - x * rate - logGamma(alpha));
+ }
+
+ @Override
+ public String toString() {
+ return this.getClass().getName() + '(' + rate + ',' + alpha + ')';
+ }
+
+ /** Returns a quick approximation of <tt>log(gamma(x))</tt>. */
+ public static double logGamma(double x) {
+
+ if (x <= 0.0 /* || x > 1.3e19 */) {
+ return -999;
+ }
+
+ double z;
+ for (z = 1.0; x < 11.0; x++) {
+ z *= x;
+ }
+
+ double r = 1.0 / (x * x);
+ double c6 = -1.9175269175269175e-03;
+ double c5 = 8.4175084175084175e-04;
+ double c4 = -5.9523809523809524e-04;
+ double c3 = 7.9365079365079365e-04;
+ double c2 = -2.7777777777777777e-03;
+ double c1 = 8.3333333333333333e-02;
+ double g = c1 + r * (c2 + r * (c3 + r * (c4 + r * (c5 + r + c6))));
+ double c0 = 9.1893853320467274e-01;
+ g = (x - 0.5) * Math.log(x) - x + c0 + g / x;
+ if (z == 1.0) {
+ return g;
+ }
+ return g - Math.log(z);
+ }
+
+}
http://git-wip-us.apache.org/repos/asf/mahout/blob/545648f6/core/src/main/java/org/apache/mahout/math/jet/random/NegativeBinomial.java
----------------------------------------------------------------------
diff --git a/core/src/main/java/org/apache/mahout/math/jet/random/NegativeBinomial.java b/core/src/main/java/org/apache/mahout/math/jet/random/NegativeBinomial.java
new file mode 100644
index 0000000..1e577eb
--- /dev/null
+++ b/core/src/main/java/org/apache/mahout/math/jet/random/NegativeBinomial.java
@@ -0,0 +1,106 @@
+/*
+ * 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.
+ */
+
+/*
+Copyright � 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
+is hereby granted without fee, provided that the above copyright notice appear in all copies and
+that both that copyright notice and this permission notice appear in supporting documentation.
+CERN makes no representations about the suitability of this software for any purpose.
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.math.jet.random;
+
+import org.apache.mahout.math.jet.math.Arithmetic;
+import org.apache.mahout.math.jet.stat.Probability;
+
+import java.util.Random;
+
+/** Mostly deprecated until unit tests are in place. Until this time, this class/interface is unsupported. */
+public final class NegativeBinomial extends AbstractDiscreteDistribution {
+
+ private final int r;
+ private final double p;
+
+ private final Gamma gamma;
+ private final Poisson poisson;
+
+ /**
+ * Constructs a Negative Binomial distribution which describes the probability of getting
+ * a particular number of negative trials (k) before getting a fixed number of positive
+ * trials (r) where each positive trial has probability (p) of being successful.
+ *
+ * @param r the required number of positive trials.
+ * @param p the probability of success.
+ * @param randomGenerator a uniform random number generator.
+ */
+ public NegativeBinomial(int r, double p, Random randomGenerator) {
+ setRandomGenerator(randomGenerator);
+ this.r = r;
+ this.p = p;
+ this.gamma = new Gamma(r, 1, randomGenerator);
+ this.poisson = new Poisson(0.0, randomGenerator);
+ }
+
+ /**
+ * Returns the cumulative distribution function.
+ */
+ public double cdf(int k) {
+ return Probability.negativeBinomial(k, r, p);
+ }
+
+ /**
+ * Returns the probability distribution function.
+ */
+ public double pdf(int k) {
+ return Arithmetic.binomial(k + r - 1, r - 1) * Math.pow(p, r) * Math.pow(1.0 - p, k);
+ }
+
+ @Override
+ public int nextInt() {
+ return nextInt(r, p);
+ }
+
+ /**
+ * Returns a sample from this distribution. The value returned will
+ * be the number of negative samples required before achieving r
+ * positive samples. Each successive sample is taken independently
+ * from a Bernouli process with probability p of success.
+ *
+ * The algorithm used is taken from J.H. Ahrens, U. Dieter (1974):
+ * Computer methods for sampling from gamma, beta, Poisson and
+ * binomial distributions, Computing 12, 223--246.
+ *
+ * This algorithm is essentially the same as described at
+ * http://en.wikipedia.org/wiki/Negative_binomial_distribution#Gamma.E2.80.93Poisson_mixture
+ * except that the notion of positive and negative outcomes is uniformly
+ * inverted. Because the inversion is complete and consistent, this
+ * definition is effectively identical to that defined on wikipedia.
+ */
+ public int nextInt(int r, double p) {
+ return this.poisson.nextInt(gamma.nextDouble(r, p / (1.0 - p)));
+ }
+
+ /**
+ * Returns a String representation of the receiver.
+ */
+ @Override
+ public String toString() {
+ return this.getClass().getName() + '(' + r + ',' + p + ')';
+ }
+
+}
http://git-wip-us.apache.org/repos/asf/mahout/blob/545648f6/core/src/main/java/org/apache/mahout/math/jet/random/Normal.java
----------------------------------------------------------------------
diff --git a/core/src/main/java/org/apache/mahout/math/jet/random/Normal.java b/core/src/main/java/org/apache/mahout/math/jet/random/Normal.java
new file mode 100644
index 0000000..7ceac22
--- /dev/null
+++ b/core/src/main/java/org/apache/mahout/math/jet/random/Normal.java
@@ -0,0 +1,110 @@
+/*
+Copyright � 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
+is hereby granted without fee, provided that the above copyright notice appear in all copies and
+that both that copyright notice and this permission notice appear in supporting documentation.
+CERN makes no representations about the suitability of this software for any purpose.
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.math.jet.random;
+
+import org.apache.mahout.math.jet.stat.Probability;
+
+import java.util.Locale;
+import java.util.Random;
+
+/**
+ * Implements a normal distribution specified mean and standard deviation.
+ */
+public class Normal extends AbstractContinousDistribution {
+
+ private double mean;
+ private double variance;
+ private double standardDeviation;
+
+ private double cache; // cache for Box-Mueller algorithm
+ private boolean cacheFilled; // Box-Mueller
+
+ private double normalizer; // performance cache
+
+ /**
+ * @param mean The mean of the resulting distribution.
+ * @param standardDeviation The standard deviation of the distribution.
+ * @param randomGenerator The random number generator to use. This can be null if you don't
+ * need to generate any numbers.
+ */
+ public Normal(double mean, double standardDeviation, Random randomGenerator) {
+ setRandomGenerator(randomGenerator);
+ setState(mean, standardDeviation);
+ }
+
+ /**
+ * Returns the cumulative distribution function.
+ */
+ @Override
+ public double cdf(double x) {
+ return Probability.normal(mean, variance, x);
+ }
+
+ /** Returns the probability density function. */
+ @Override
+ public double pdf(double x) {
+ double diff = x - mean;
+ return normalizer * Math.exp(-(diff * diff) / (2.0 * variance));
+ }
+
+ /**
+ * Returns a random number from the distribution.
+ */
+ @Override
+ public double nextDouble() {
+ // Uses polar Box-Muller transformation.
+ if (cacheFilled) {
+ cacheFilled = false;
+ return cache;
+ }
+
+ double x;
+ double y;
+ double r;
+ do {
+ x = 2.0 * randomDouble() - 1.0;
+ y = 2.0 * randomDouble() - 1.0;
+ r = x * x + y * y;
+ } while (r >= 1.0);
+
+ double z = Math.sqrt(-2.0 * Math.log(r) / r);
+ cache = this.mean + this.standardDeviation * x * z;
+ cacheFilled = true;
+ return this.mean + this.standardDeviation * y * z;
+ }
+
+ /** Sets the uniform random generator internally used. */
+ @Override
+ public final void setRandomGenerator(Random randomGenerator) {
+ super.setRandomGenerator(randomGenerator);
+ this.cacheFilled = false;
+ }
+
+ /**
+ * Sets the mean and variance.
+ * @param mean The new value for the mean.
+ * @param standardDeviation The new value for the standard deviation.
+ */
+ public final void setState(double mean, double standardDeviation) {
+ if (mean != this.mean || standardDeviation != this.standardDeviation) {
+ this.mean = mean;
+ this.standardDeviation = standardDeviation;
+ this.variance = standardDeviation * standardDeviation;
+ this.cacheFilled = false;
+
+ this.normalizer = 1.0 / Math.sqrt(2.0 * Math.PI * variance);
+ }
+ }
+
+ /** Returns a String representation of the receiver. */
+ @Override
+ public String toString() {
+ return String.format(Locale.ENGLISH, "%s(m=%f, sd=%f)", this.getClass().getName(), mean, standardDeviation);
+ }
+}
http://git-wip-us.apache.org/repos/asf/mahout/blob/545648f6/core/src/main/java/org/apache/mahout/math/jet/random/Poisson.java
----------------------------------------------------------------------
diff --git a/core/src/main/java/org/apache/mahout/math/jet/random/Poisson.java b/core/src/main/java/org/apache/mahout/math/jet/random/Poisson.java
new file mode 100644
index 0000000..497691e
--- /dev/null
+++ b/core/src/main/java/org/apache/mahout/math/jet/random/Poisson.java
@@ -0,0 +1,296 @@
+/*
+Copyright � 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
+is hereby granted without fee, provided that the above copyright notice appear in all copies and
+that both that copyright notice and this permission notice appear in supporting documentation.
+CERN makes no representations about the suitability of this software for any purpose.
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.math.jet.random;
+
+import org.apache.mahout.math.jet.math.Arithmetic;
+
+import java.util.Random;
+
+/** Partially deprecated until unit tests are in place. Until this time, this class/interface is unsupported. */
+public final class Poisson extends AbstractDiscreteDistribution {
+
+ private final double mean;
+
+ // precomputed and cached values (for performance only)
+ // cache for < SWITCH_MEAN
+ private double myOld = -1.0;
+ private double p;
+ private double q;
+ private double p0;
+ private final double[] pp = new double[36];
+ private int llll;
+
+ // cache for >= SWITCH_MEAN
+ private double myLast = -1.0;
+ private double ll;
+ private int k2;
+ private int k4;
+ private int k1;
+ private int k5;
+ private double dl;
+ private double dr;
+ private double r1;
+ private double r2;
+ private double r4;
+ private double r5;
+ private double lr;
+ private double lMy;
+ private double cPm;
+ private double f1;
+ private double f2;
+ private double f4;
+ private double f5;
+ private double p1;
+ private double p2;
+ private double p3;
+ private double p4;
+ private double p5;
+ private double p6;
+
+ // cache for both;
+
+
+ private static final double MEAN_MAX = Integer.MAX_VALUE;
+ // for all means larger than that, we don't try to compute a poisson deviation, but return the mean.
+ private static final double SWITCH_MEAN = 10.0; // switch from method A to method B
+
+
+ /** Constructs a poisson distribution. Example: mean=1.0. */
+ public Poisson(double mean, Random randomGenerator) {
+ setRandomGenerator(randomGenerator);
+ this.mean = mean;
+ }
+
+ private static double f(int k, double lNu, double cPm) {
+ return Math.exp(k * lNu - Arithmetic.logFactorial(k) - cPm);
+ }
+
+ @Override
+ public int nextInt() {
+ return nextInt(mean);
+ }
+
+ /** Returns a random number from the distribution; bypasses the internal state. */
+ public int nextInt(double theMean) {
+ /******************************************************************
+ * *
+ * Poisson Distribution - Patchwork Rejection/Inversion *
+ * *
+ ******************************************************************
+ * *
+ * For parameter my < 10 Tabulated Inversion is applied. *
+ * For my >= 10 Patchwork Rejection is employed: *
+ * The area below the histogram function f(x) is rearranged in *
+ * its body by certain point reflections. Within a large center *
+ * interval variates are sampled efficiently by rejection from *
+ * uniform hats. Rectangular immediate acceptance regions speed *
+ * up the generation. The remaining tails are covered by *
+ * exponential functions. *
+ * *
+ *****************************************************************/
+ Random gen = getRandomGenerator();
+
+ //double t, g, my_k;
+
+ //double gx, gy, px, py, e, x, xx, delta, v;
+ //int sign;
+
+ //static double p,q,p0,pp[36];
+ //static long ll,m;
+
+ int m;
+ if (theMean < SWITCH_MEAN) { // CASE B: Inversion- start new table and calculate p0
+ if (theMean != myOld) {
+ myOld = theMean;
+ llll = 0;
+ p = Math.exp(-theMean);
+ q = p;
+ p0 = p;
+ //for (k=pp.length; --k >=0;) pp[k] = 0;
+ }
+ m = theMean > 1.0 ? (int) theMean : 1;
+ while (true) {
+ double u = gen.nextDouble();
+ int k = 0;
+ if (u <= p0) {
+ return k;
+ }
+ if (llll != 0) { // Step T. Table comparison
+ int i = u > 0.458 ? Math.min(llll, m) : 1;
+ for (k = i; k <= llll; k++) {
+ if (u <= pp[k]) {
+ return k;
+ }
+ }
+ if (llll == 35) {
+ continue;
+ }
+ }
+ for (k = llll + 1; k <= 35; k++) { // Step C. Creation of new prob.
+ p *= theMean / k;
+ q += p;
+ pp[k] = q;
+ if (u <= q) {
+ llll = k;
+ return k;
+ }
+ }
+ llll = 35;
+ }
+ // end my < SWITCH_MEAN
+ } else if (theMean < MEAN_MAX) { // CASE A: acceptance complement
+ //static double my_last = -1.0;
+ //static long int m, k2, k4, k1, k5;
+ //static double dl, dr, r1, r2, r4, r5, ll, lr, l_my, c_pm,
+ // f1, f2, f4, f5, p1, p2, p3, p4, p5, p6;
+
+ m = (int) theMean;
+ if (theMean != myLast) { // set-up
+ myLast = theMean;
+
+ // approximate deviation of reflection points k2, k4 from my - 1/2
+ double Ds = Math.sqrt(theMean + 0.25);
+
+ // mode m, reflection points k2 and k4, and points k1 and k5, which
+ // delimit the centre region of h(x)
+ k2 = (int) Math.ceil(theMean - 0.5 - Ds);
+ k4 = (int) (theMean - 0.5 + Ds);
+ k1 = k2 + k2 - m + 1;
+ k5 = k4 + k4 - m;
+
+ // range width of the critical left and right centre region
+ dl = k2 - k1;
+ dr = k5 - k4;
+
+ // recurrence constants r(k) = p(k)/p(k-1) at k = k1, k2, k4+1, k5+1
+ r1 = theMean / k1;
+ r2 = theMean / k2;
+ r4 = theMean / (k4 + 1);
+ r5 = theMean / (k5 + 1);
+
+ // reciprocal values of the scale parameters of expon. tail envelopes
+ ll = Math.log(r1); // expon. tail left
+ lr = -Math.log(r5); // expon. tail right
+
+ // Poisson constants, necessary for computing function values f(k)
+ lMy = Math.log(theMean);
+ cPm = m * lMy - Arithmetic.logFactorial(m);
+
+ // function values f(k) = p(k)/p(m) at k = k2, k4, k1, k5
+ f2 = f(k2, lMy, cPm);
+ f4 = f(k4, lMy, cPm);
+ f1 = f(k1, lMy, cPm);
+ f5 = f(k5, lMy, cPm);
+
+ // area of the two centre and the two exponential tail regions
+ // area of the two immediate acceptance regions between k2, k4
+ p1 = f2 * (dl + 1.0); // immed. left
+ p2 = f2 * dl + p1; // centre left
+ p3 = f4 * (dr + 1.0) + p2; // immed. right
+ p4 = f4 * dr + p3; // centre right
+ p5 = f1 / ll + p4; // expon. tail left
+ p6 = f5 / lr + p5; // expon. tail right
+ } // end set-up
+
+ while (true) {
+ // generate uniform number U -- U(0, p6)
+ // case distinction corresponding to U
+ double W;
+ double V;
+ double U;
+ int Y;
+ int X;
+ int Dk;
+ if ((U = gen.nextDouble() * p6) < p2) { // centre left
+
+ // immediate acceptance region R2 = [k2, m) *[0, f2), X = k2, ... m -1
+ if ((V = U - p1) < 0.0) {
+ return k2 + (int) (U / f2);
+ }
+ // immediate acceptance region R1 = [k1, k2)*[0, f1), X = k1, ... k2-1
+ if ((W = V / dl) < f1) {
+ return k1 + (int) (V / f1);
+ }
+
+ // computation of candidate X < k2, and its counterpart Y > k2
+ // either squeeze-acceptance of X or acceptance-rejection of Y
+ Dk = gen.nextInt((int) dl) + 1;
+ if (W <= f2 - Dk * (f2 - f2 / r2)) { // quick accept of
+ return k2 - Dk; // X = k2 - Dk
+ }
+ if ((V = f2 + f2 - W) < 1.0) { // quick reject of Y
+ Y = k2 + Dk;
+ if (V <= f2 + Dk * (1.0 - f2) / (dl + 1.0)) { // quick accept of
+ return Y; // Y = k2 + Dk
+ }
+ if (V <= f(Y, lMy, cPm)) {
+ return Y;
+ } // final accept of Y
+ }
+ X = k2 - Dk;
+ } else if (U < p4) { // centre right
+ // immediate acceptance region R3 = [m, k4+1)*[0, f4), X = m, ... k4
+ if ((V = U - p3) < 0.0) {
+ return k4 - (int) ((U - p2) / f4);
+ }
+ // immediate acceptance region R4 = [k4+1, k5+1)*[0, f5)
+ if ((W = V / dr) < f5) {
+ return k5 - (int) (V / f5);
+ }
+
+ // computation of candidate X > k4, and its counterpart Y < k4
+ // either squeeze-acceptance of X or acceptance-rejection of Y
+ Dk = gen.nextInt((int) dr) + 1;
+ if (W <= f4 - Dk * (f4 - f4 * r4)) { // quick accept of
+ return k4 + Dk; // X = k4 + Dk
+ }
+ if ((V = f4 + f4 - W) < 1.0) { // quick reject of Y
+ Y = k4 - Dk;
+ if (V <= f4 + Dk * (1.0 - f4) / dr) { // quick accept of
+ return Y; // Y = k4 - Dk
+ }
+ if (V <= f(Y, lMy, cPm)) {
+ return Y;
+ } // final accept of Y
+ }
+ X = k4 + Dk;
+ } else {
+ W = gen.nextDouble();
+ if (U < p5) { // expon. tail left
+ Dk = (int) (1.0 - Math.log(W) / ll);
+ if ((X = k1 - Dk) < 0) {
+ continue;
+ } // 0 <= X <= k1 - 1
+ W *= (U - p4) * ll; // W -- U(0, h(x))
+ if (W <= f1 - Dk * (f1 - f1 / r1)) {
+ return X;
+ } // quick accept of X
+ } else { // expon. tail right
+ Dk = (int) (1.0 - Math.log(W) / lr);
+ X = k5 + Dk; // X >= k5 + 1
+ W *= (U - p5) * lr; // W -- U(0, h(x))
+ if (W <= f5 - Dk * (f5 - f5 * r5)) {
+ return X;
+ } // quick accept of X
+ }
+ }
+
+ // acceptance-rejection test of candidate X from the original area
+ // test, whether W <= f(k), with W = U*h(x) and U -- U(0, 1)
+ // log f(X) = (X - m)*log(my) - log X! + log m!
+ if (Math.log(W) <= X * lMy - Arithmetic.logFactorial(X) - cPm) {
+ return X;
+ }
+ }
+ } else { // mean is too large
+ return (int) theMean;
+ }
+ }
+
+}
http://git-wip-us.apache.org/repos/asf/mahout/blob/545648f6/core/src/main/java/org/apache/mahout/math/jet/random/Uniform.java
----------------------------------------------------------------------
diff --git a/core/src/main/java/org/apache/mahout/math/jet/random/Uniform.java b/core/src/main/java/org/apache/mahout/math/jet/random/Uniform.java
new file mode 100644
index 0000000..32c8b90
--- /dev/null
+++ b/core/src/main/java/org/apache/mahout/math/jet/random/Uniform.java
@@ -0,0 +1,164 @@
+/*
+Copyright � 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
+is hereby granted without fee, provided that the above copyright notice appear in all copies and
+that both that copyright notice and this permission notice appear in supporting documentation.
+CERN makes no representations about the suitability of this software for any purpose.
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.math.jet.random;
+
+import org.apache.mahout.common.RandomUtils;
+
+import java.util.Random;
+
+public class Uniform extends AbstractContinousDistribution {
+
+ private double min;
+ private double max;
+
+ /**
+ * Constructs a uniform distribution with the given minimum and maximum, using a {@link
+ * org.apache.mahout.math.jet.random.engine.MersenneTwister} seeded with the given seed.
+ */
+ public Uniform(double min, double max, int seed) {
+ this(min, max, RandomUtils.getRandom(seed));
+ }
+
+ /** Constructs a uniform distribution with the given minimum and maximum. */
+ public Uniform(double min, double max, Random randomGenerator) {
+ setRandomGenerator(randomGenerator);
+ setState(min, max);
+ }
+
+ /** Constructs a uniform distribution with <tt>min=0.0</tt> and <tt>max=1.0</tt>. */
+ public Uniform(Random randomGenerator) {
+ this(0, 1, randomGenerator);
+ }
+
+ /** Returns the cumulative distribution function (assuming a continous uniform distribution). */
+ @Override
+ public double cdf(double x) {
+ if (x <= min) {
+ return 0.0;
+ }
+ if (x >= max) {
+ return 1.0;
+ }
+ return (x - min) / (max - min);
+ }
+
+ /** Returns a uniformly distributed random <tt>boolean</tt>. */
+ public boolean nextBoolean() {
+ return randomDouble() > 0.5;
+ }
+
+ /**
+ * Returns a uniformly distributed random number in the open interval <tt>(min,max)</tt> (excluding <tt>min</tt> and
+ * <tt>max</tt>).
+ */
+ @Override
+ public double nextDouble() {
+ return min + (max - min) * randomDouble();
+ }
+
+ /**
+ * Returns a uniformly distributed random number in the open interval <tt>(from,to)</tt> (excluding <tt>from</tt> and
+ * <tt>to</tt>). Pre conditions: <tt>from <= to</tt>.
+ */
+ public double nextDoubleFromTo(double from, double to) {
+ return from + (to - from) * randomDouble();
+ }
+
+ /**
+ * Returns a uniformly distributed random number in the open interval <tt>(from,to)</tt> (excluding <tt>from</tt> and
+ * <tt>to</tt>). Pre conditions: <tt>from <= to</tt>.
+ */
+ public float nextFloatFromTo(float from, float to) {
+ return (float) nextDoubleFromTo(from, to);
+ }
+
+ /**
+ * Returns a uniformly distributed random number in the closed interval
+ * <tt>[from,to]</tt> (including <tt>from</tt>
+ * and <tt>to</tt>). Pre conditions: <tt>from <= to</tt>.
+ */
+ public int nextIntFromTo(int from, int to) {
+ return (int) (from + (long) ((1L + to - from) * randomDouble()));
+ }
+
+ /**
+ * Returns a uniformly distributed random number in the closed interval <tt>[from,to]</tt> (including <tt>from</tt>
+ * and <tt>to</tt>). Pre conditions: <tt>from <= to</tt>.
+ */
+ public long nextLongFromTo(long from, long to) {
+ /* Doing the thing turns out to be more tricky than expected.
+ avoids overflows and underflows.
+ treats cases like from=-1, to=1 and the like right.
+ the following code would NOT solve the problem: return (long) (Doubles.randomFromTo(from,to));
+
+ rounding avoids the unsymmetric behaviour of casts from double to long: (long) -0.7 = 0, (long) 0.7 = 0.
+ checking for overflows and underflows is also necessary.
+ */
+
+ // first the most likely and also the fastest case.
+ if (from >= 0 && to < Long.MAX_VALUE) {
+ return from + (long) nextDoubleFromTo(0.0, to - from + 1);
+ }
+
+ // would we get a numeric overflow?
+ // if not, we can still handle the case rather efficient.
+ double diff = (double) to - (double) from + 1.0;
+ if (diff <= Long.MAX_VALUE) {
+ return from + (long) nextDoubleFromTo(0.0, diff);
+ }
+
+ // now the pathologic boundary cases.
+ // they are handled rather slow.
+ long random;
+ if (from == Long.MIN_VALUE) {
+ if (to == Long.MAX_VALUE) {
+ //return Math.round(nextDoubleFromTo(from,to));
+ int i1 = nextIntFromTo(Integer.MIN_VALUE, Integer.MAX_VALUE);
+ int i2 = nextIntFromTo(Integer.MIN_VALUE, Integer.MAX_VALUE);
+ return ((i1 & 0xFFFFFFFFL) << 32) | (i2 & 0xFFFFFFFFL);
+ }
+ random = Math.round(nextDoubleFromTo(Long.MIN_VALUE, to + 1));
+ if (random > to) {
+ random = Long.MIN_VALUE;
+ }
+ } else {
+ random = Math.round(nextDoubleFromTo(from - 1, to));
+ if (random < from) {
+ random = to;
+ }
+ }
+ return random;
+ }
+
+ /** Returns the probability distribution function (assuming a continous uniform distribution). */
+ @Override
+ public double pdf(double x) {
+ if (x <= min || x >= max) {
+ return 0.0;
+ }
+ return 1.0 / (max - min);
+ }
+
+ /** Sets the internal state. */
+ public void setState(double min, double max) {
+ if (max < min) {
+ setState(max, min);
+ return;
+ }
+ this.min = min;
+ this.max = max;
+ }
+
+
+ /** Returns a String representation of the receiver. */
+ @Override
+ public String toString() {
+ return this.getClass().getName() + '(' + min + ',' + max + ')';
+ }
+}
http://git-wip-us.apache.org/repos/asf/mahout/blob/545648f6/core/src/main/java/org/apache/mahout/math/jet/random/engine/MersenneTwister.java
----------------------------------------------------------------------
diff --git a/core/src/main/java/org/apache/mahout/math/jet/random/engine/MersenneTwister.java b/core/src/main/java/org/apache/mahout/math/jet/random/engine/MersenneTwister.java
new file mode 100644
index 0000000..8bca895
--- /dev/null
+++ b/core/src/main/java/org/apache/mahout/math/jet/random/engine/MersenneTwister.java
@@ -0,0 +1,275 @@
+/**
+ * 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.
+ */
+
+/*
+Copyright 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
+is hereby granted without fee, provided that the above copyright notice appear in all copies and
+that both that copyright notice and this permission notice appear in supporting documentation.
+CERN makes no representations about the suitability of this software for any purpose.
+It is provided "as is" without expressed or implied warranty.
+*/
+
+package org.apache.mahout.math.jet.random.engine;
+
+import java.util.Date;
+/**
+ MersenneTwister (MT19937) is one of the strongest uniform pseudo-random number generators
+ known so far; at the same time it is quick.
+ Produces uniformly distributed <tt>int</tt>'s and <tt>long</tt>'s in the closed intervals
+ <tt>[Integer.MIN_VALUE,Integer.MAX_VALUE]</tt> and <tt>[Long.MIN_VALUE,Long.MAX_VALUE]</tt>, respectively,
+ as well as <tt>float</tt>'s and <tt>double</tt>'s in the open unit intervals <tt>(0.0f,1.0f)</tt>
+ and <tt>(0.0,1.0)</tt>, respectively.
+ The seed can be any 32-bit integer except <tt>0</tt>. Shawn J. Cokus commented that perhaps the
+ seed should preferably be odd.
+ <p>
+ <b>Quality:</b> MersenneTwister is designed to pass the k-distribution test. It has an
+ astronomically large period of 2<sup>19937</sup>-1 (=10<sup>6001</sup>) and 623-dimensional
+ equidistribution up to 32-bit accuracy.
+ It passes many stringent statistical tests, including the
+ <A HREF="http://stat.fsu.edu/~geo/diehard.html">diehard</A> test of G. Marsaglia
+ and the load test of P. Hellekalek and S. Wegenkittl.
+ <p>
+ <b>Performance:</b> Its speed is comparable to other modern generators (in particular,
+ as fast as <tt>java.util.Random.nextFloat()</tt>).
+ 2.5 million calls to <tt>raw()</tt> per second (Pentium Pro 200 Mhz, JDK 1.2, NT).
+ Be aware, however, that there is a non-negligible amount of overhead required to initialize the data
+ structures used by a MersenneTwister. Code like
+ {@code
+ double sum = 0.0;
+ for (int i=0; i<100000; ++i) {
+ RandomElement twister = new MersenneTwister(new Date());
+ sum += twister.raw();
+ }
+ }
+ will be wildly inefficient. Consider using
+ {@code
+ double sum = 0.0;
+ RandomElement twister = new MersenneTwister(new Date());
+ for (int i=0; i<100000; ++i) {
+ sum += twister.raw();
+ }
+ }
+ instead. This allows the cost of constructing the MersenneTwister object
+ to be borne only once, rather than once for each iteration in the loop.
+ <p>
+ <b>Implementation:</b> After M. Matsumoto and T. Nishimura,
+ "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform Pseudo-Random Number Generator",
+ ACM Transactions on Modeling and Computer Simulation,
+ Vol. 8, No. 1, January 1998, pp 3--30.
+ <dl>
+ <dt>More info on <a HREF="http://www.math.keio.ac.jp/~matumoto/eindex.html"> Masumoto's homepage</a>.</dt>
+ <dt>More info on <a HREF="http://www.ncsa.uiuc.edu/Apps/CMP/RNG/www-rng.html"> Pseudo-random number
+ generators is on the Web</a>.</dt>
+ <dt>Yet <a HREF="http://nhse.npac.syr.edu/random"> some more info</a>.</dt>
+ <p>
+ The correctness of this implementation has been verified against the published output sequence
+ <a href="http://www.math.keio.ac.jp/~nisimura/random/real2/mt19937-2.out">mt19937-2.out</a> of the C-implementation
+ <a href="http://www.math.keio.ac.jp/~nisimura/random/real2/mt19937-2.c">mt19937-2.c</a>.
+ (Call <tt>test(1000)</tt> to print the sequence).
+ <dt>
+ Note that this implementation is <b>not synchronized</b>.</dt>
+ </dl>
+ <p>
+ <b>Details:</b> MersenneTwister is designed with consideration of the flaws of various existing generators in mind.
+ It is an improved version of TT800, a very successful generator.
+ MersenneTwister is based on linear recurrences modulo 2.
+ Such generators are very fast, have extremely long periods, and appear quite robust.
+ MersenneTwister produces 32-bit numbers, and every <tt>k</tt>-dimensional vector of such
+ numbers appears the same number of times as <tt>k</tt> successive values over the
+ period length, for each <tt>k <= 623</tt> (except for the zero vector, which appears one time less).
+ If one looks at only the first <tt>n <= 16</tt> bits of each number, then the property holds
+ for even larger <tt>k</tt>, as shown in the following table (taken from the publication cited above):
+ <table width="75%" border="1" cellspacing="0" cellpadding="0" summary="property table">
+ <tr>
+ <td width="2%" align="center"> <div>n</div> </td>
+ <td width="6%" align="center"> <div>1</div> </td>
+ <td width="5%" align="center"> <div>2</div> </td>
+ <td width="5%" align="center"> <div>3</div> </td>
+ <td width="5%" align="center"> <div>4</div> </td>
+ <td width="5%" align="center"> <div>5</div> </td>
+ <td width="5%" align="center"> <div>6</div> </td>
+ <td width="5%" align="center"> <div>7</div> </td>
+ <td width="5%" align="center"> <div>8</div> </td>
+ <td width="5%" align="center"> <div>9</div> </td>
+ <td width="5%" align="center"> <div>10</div> </td>
+ <td width="5%" align="center"> <div>11</div> </td>
+ <td width="10%" align="center"> <div>12 .. 16</div> </td>
+ <td width="10%" align="center"> <div>17 .. 32</div> </td>
+ </tr>
+ <tr>
+ <td width="2%" align="center"> <div>k</div> </td>
+ <td width="6%" align="center"> <div>19937</div> </td>
+ <td width="5%" align="center"> <div>9968</div> </td>
+ <td width="5%" align="center"> <div>6240</div> </td>
+ <td width="5%" align="center"> <div>4984</div> </td>
+ <td width="5%" align="center"> <div>3738</div> </td>
+ <td width="5%" align="center"> <div>3115</div> </td>
+ <td width="5%" align="center"> <div>2493</div> </td>
+ <td width="5%" align="center"> <div>2492</div> </td>
+ <td width="5%" align="center"> <div>1869</div> </td>
+ <td width="5%" align="center"> <div>1869</div> </td>
+ <td width="5%" align="center"> <div>1248</div> </td>
+ <td width="10%" align="center"> <div>1246</div> </td>
+ <td width="10%" align="center"> <div>623</div> </td>
+ </tr>
+ </table>
+ <p>
+ MersenneTwister generates random numbers in batches of 624 numbers at a time, so
+ the caching and pipelining of modern systems is exploited.
+ The generator is implemented to generate the output by using the fastest arithmetic
+ operations only: 32-bit additions and bit operations (no division, no multiplication, no mod).
+ These operations generate sequences of 32 random bits (<tt>int</tt>'s).
+ <tt>long</tt>'s are formed by concatenating two 32 bit <tt>int</tt>'s.
+ <tt>float</tt>'s are formed by dividing the interval <tt>[0.0,1.0]</tt> into 2<sup>32</sup>
+ sub intervals, then randomly choosing one subinterval.
+ <tt>double</tt>'s are formed by dividing the interval <tt>[0.0,1.0]</tt> into 2<sup>64</sup>
+ sub intervals, then randomly choosing one subinterval.
+ <p>
+ @author wolfgang.hoschek@cern.ch
+ @version 1.0, 09/24/99
+ @see java.util.Random
+ */
+public final class MersenneTwister extends RandomEngine {
+
+ /* Period parameters */
+ private static final int N = 624;
+ private static final int M = 397;
+ private static final int MATRIX_A = 0x9908b0df; /* constant vector a */
+ private static final int UPPER_MASK = 0x80000000; /* most significant w-r bits */
+ private static final int LOWER_MASK = 0x7fffffff; /* least significant r bits */
+
+ /* for tempering */
+ private static final int TEMPERING_MASK_B = 0x9d2c5680;
+ private static final int TEMPERING_MASK_C = 0xefc60000;
+
+ private static final int MAG0 = 0x0;
+ private static final int MAG1 = MATRIX_A;
+ //private static final int[] mag01=new int[] {0x0, MATRIX_A};
+ /* mag01[x] = x * MATRIX_A for x=0,1 */
+
+ private static final int DEFAULT_SEED = 4357;
+
+ private int mti;
+ private final int[] mt = new int[N]; /* set initial seeds: N = 624 words */
+
+ /**
+ * Constructs and returns a random number generator with a default seed, which is a <b>constant</b>. Thus using this
+ * constructor will yield generators that always produce exactly the same sequence. This method is mainly intended to
+ * ease testing and debugging.
+ */
+ public MersenneTwister() {
+ this(DEFAULT_SEED);
+ }
+
+ /** Constructs and returns a random number generator with the given seed.
+ * @param seed A number that is used to initialize the internal state of the generator.
+ */
+ public MersenneTwister(int seed) {
+ setSeed(seed);
+ }
+
+ /**
+ * Constructs and returns a random number generator seeded with the given date.
+ *
+ * @param d typically <tt>new Date()</tt>
+ */
+ public MersenneTwister(Date d) {
+ this((int) d.getTime());
+ }
+
+ /** Generates N words at one time */
+ void nextBlock() {
+ int y;
+ int kk;
+
+ for (kk = 0; kk < N - M; kk++) {
+ y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
+ mt[kk] = mt[kk + M] ^ (y >>> 1) ^ ((y & 0x1) == 0 ? MAG0 : MAG1);
+ }
+ for (; kk < N - 1; kk++) {
+ y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
+ mt[kk] = mt[kk + (M - N)] ^ (y >>> 1) ^ ((y & 0x1) == 0 ? MAG0 : MAG1);
+ }
+ y = (mt[N - 1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
+ mt[N - 1] = mt[M - 1] ^ (y >>> 1) ^ ((y & 0x1) == 0 ? MAG0 : MAG1);
+
+ this.mti = 0;
+ }
+
+ /**
+ * Returns a 32 bit uniformly distributed random number in the closed interval
+ * <tt>[Integer.MIN_VALUE,Integer.MAX_VALUE]</tt>
+ * (including <tt>Integer.MIN_VALUE</tt> and <tt>Integer.MAX_VALUE</tt>).
+ */
+ @Override
+ public int nextInt() {
+ /* Each single bit including the sign bit will be random */
+ if (mti == N) {
+ nextBlock();
+ } // generate N ints at one time
+
+ int y = mt[mti++];
+ y ^= y >>> 11; // y ^= TEMPERING_SHIFT_U(y );
+ y ^= (y << 7) & TEMPERING_MASK_B; // y ^= TEMPERING_SHIFT_S(y) & TEMPERING_MASK_B;
+ y ^= (y << 15) & TEMPERING_MASK_C; // y ^= TEMPERING_SHIFT_T(y) & TEMPERING_MASK_C;
+ // y &= 0xffffffff; //you may delete this line if word size = 32
+ y ^= y >>> 18; // y ^= TEMPERING_SHIFT_L(y);
+
+ return y;
+ }
+
+ /** Sets the receiver's seed. This method resets the receiver's entire internal state.
+ * @param seed An integer that is used to reset the internal state of the generator */
+ void setSeed(int seed) {
+ mt[0] = seed;
+ for (int i = 1; i < N; i++) {
+ mt[i] = 1812433253 * (mt[i - 1] ^ (mt[i - 1] >> 30)) + i;
+ /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
+ /* In the previous versions, MSBs of the seed affect */
+ /* only MSBs of the array mt[]. */
+ /* 2002/01/09 modified by Makoto Matsumoto */
+ //mt[i] &= 0xffffffff;
+ /* for >32 bit machines */
+ }
+ //log.info("init done");
+ mti = N;
+ }
+
+ /**
+ * Sets the receiver's seed in a fashion compatible with the
+ * reference C implementation. See
+ * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/VERSIONS/C-LANG/980409/mt19937int.c
+ *
+ * This method isn't as good as the default method due to poor distribution of the
+ * resulting states.
+ *
+ * @param seed An integer that is used to reset the internal state in the same way as
+ * done in the 1999 reference implementation. Should only be used for testing, not
+ * actual coding.
+ */
+ void setReferenceSeed(int seed) {
+ for (int i = 0; i < N; i++) {
+ mt[i] = seed & 0xffff0000;
+ seed = 69069 * seed + 1;
+ mt[i] |= (seed & 0xffff0000) >>> 16;
+ seed = 69069 * seed + 1;
+ }
+ //log.info("init done");
+ mti = N;
+ }
+}
http://git-wip-us.apache.org/repos/asf/mahout/blob/545648f6/core/src/main/java/org/apache/mahout/math/jet/random/engine/RandomEngine.java
----------------------------------------------------------------------
diff --git a/core/src/main/java/org/apache/mahout/math/jet/random/engine/RandomEngine.java b/core/src/main/java/org/apache/mahout/math/jet/random/engine/RandomEngine.java
new file mode 100644
index 0000000..f832b1d
--- /dev/null
+++ b/core/src/main/java/org/apache/mahout/math/jet/random/engine/RandomEngine.java
@@ -0,0 +1,169 @@
+/**
+ * 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.
+ */
+/*
+Copyright � 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
+is hereby granted without fee, provided that the above copyright notice appear in all copies and
+that both that copyright notice and this permission notice appear in supporting documentation.
+CERN makes no representations about the suitability of this software for any purpose.
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.math.jet.random.engine;
+
+import org.apache.mahout.math.function.DoubleFunction;
+import org.apache.mahout.math.function.IntFunction;
+
+/**
+ * Abstract base class for uniform pseudo-random number generating engines.
+ * <p>
+ * Most probability distributions are obtained by using a <b>uniform</b> pseudo-random number generation engine
+ * followed by a transformation to the desired distribution.
+ * Thus, subclasses of this class are at the core of computational statistics, simulations, Monte Carlo methods, etc.
+ * <p>
+ * Subclasses produce uniformly distributed <tt>int</tt>'s and <tt>long</tt>'s in the closed intervals
+ * <tt>[Integer.MIN_VALUE,Integer.MAX_VALUE]</tt> and <tt>[Long.MIN_VALUE,Long.MAX_VALUE]</tt>, respectively,
+ * as well as <tt>float</tt>'s and <tt>double</tt>'s in the open unit intervals <tt>(0.0f,1.0f)</tt> and
+ * <tt>(0.0,1.0)</tt>, respectively.
+ * <p>
+ * Subclasses need to override one single method only: <tt>nextInt()</tt>.
+ * All other methods generating different data types or ranges are usually layered upon <tt>nextInt()</tt>.
+ * <tt>long</tt>'s are formed by concatenating two 32 bit <tt>int</tt>'s.
+ * <tt>float</tt>'s are formed by dividing the interval <tt>[0.0f,1.0f]</tt> into 2<sup>32</sup> sub intervals,
+ * then randomly choosing one subinterval.
+ * <tt>double</tt>'s are formed by dividing the interval <tt>[0.0,1.0]</tt> into 2<sup>64</sup> sub intervals,
+ * then randomly choosing one subinterval.
+ * <p>
+ * Note that this implementation is <b>not synchronized</b>.
+ *
+ * @see MersenneTwister
+ * @see java.util.Random
+ */
+public abstract class RandomEngine extends DoubleFunction implements IntFunction {
+
+ /**
+ * Equivalent to <tt>raw()</tt>. This has the effect that random engines can now be used as function objects,
+ * returning a random number upon function evaluation.
+ */
+ @Override
+ public double apply(double dummy) {
+ return raw();
+ }
+
+ /**
+ * Equivalent to <tt>nextInt()</tt>. This has the effect that random engines can now be used as function objects,
+ * returning a random number upon function evaluation.
+ */
+ @Override
+ public int apply(int dummy) {
+ return nextInt();
+ }
+
+ /**
+ * @return a 64 bit uniformly distributed random number in the open unit interval {@code (0.0,1.0)} (excluding
+ * 0.0 and 1.0).
+ */
+ public double nextDouble() {
+ double nextDouble;
+
+ do {
+ // -9.223372036854776E18 == (double) Long.MIN_VALUE
+ // 5.421010862427522E-20 == 1 / Math.pow(2,64) == 1 / ((double) Long.MAX_VALUE - (double) Long.MIN_VALUE);
+ nextDouble = (nextLong() - -9.223372036854776E18) * 5.421010862427522E-20;
+ }
+ // catch loss of precision of long --> double conversion
+ while (!(nextDouble > 0.0 && nextDouble < 1.0));
+
+ // --> in (0.0,1.0)
+ return nextDouble;
+
+ /*
+ nextLong == Long.MAX_VALUE --> 1.0
+ nextLong == Long.MIN_VALUE --> 0.0
+ nextLong == Long.MAX_VALUE-1 --> 1.0
+ nextLong == Long.MAX_VALUE-100000L --> 0.9999999999999946
+ nextLong == Long.MIN_VALUE+1 --> 0.0
+ nextLong == Long.MIN_VALUE-100000L --> 0.9999999999999946
+ nextLong == 1L --> 0.5
+ nextLong == -1L --> 0.5
+ nextLong == 2L --> 0.5
+ nextLong == -2L --> 0.5
+ nextLong == 2L+100000L --> 0.5000000000000054
+ nextLong == -2L-100000L --> 0.49999999999999456
+ */
+ }
+
+ /**
+ * @return a 32 bit uniformly distributed random number in the open unit interval {@code (0.0f, 1.0f)} (excluding
+ * 0.0f and 1.0f).
+ */
+ public float nextFloat() {
+ // catch loss of precision of double --> float conversion which could result in a value == 1.0F
+ float nextFloat;
+ do {
+ nextFloat = (float) raw();
+ }
+ while (nextFloat >= 1.0f);
+
+ // --> in [0.0f,1.0f)
+ return nextFloat;
+ }
+
+ /**
+ * @return a 32 bit uniformly distributed random number in the closed interval
+ * <tt>[Integer.MIN_VALUE,Integer.MAX_VALUE]</tt>
+ * (including <tt>Integer.MIN_VALUE</tt> and <tt>Integer.MAX_VALUE</tt>);
+ */
+ public abstract int nextInt();
+
+ /**
+ * @return a 64 bit uniformly distributed random number in the closed interval
+ * <tt>[Long.MIN_VALUE,Long.MAX_VALUE]</tt>
+ * (including <tt>Long.MIN_VALUE</tt> and <tt>Long.MAX_VALUE</tt>).
+ */
+ public long nextLong() {
+ // concatenate two 32-bit strings into one 64-bit string
+ return ((nextInt() & 0xFFFFFFFFL) << 32) | (nextInt() & 0xFFFFFFFFL);
+ }
+
+ /**
+ * @return a 32 bit uniformly distributed random number in the open unit interval {@code (0.0, 1.0)} (excluding
+ * 0.0 and 1.0).
+ */
+ public double raw() {
+ int nextInt;
+ do { // accept anything but zero
+ nextInt = nextInt(); // in [Integer.MIN_VALUE,Integer.MAX_VALUE]-interval
+ } while (nextInt == 0);
+
+ // transform to (0.0,1.0)-interval
+ // 2.3283064365386963E-10 == 1.0 / Math.pow(2,32)
+ return (nextInt & 0xFFFFFFFFL) * 2.3283064365386963E-10;
+
+ /*
+ nextInt == Integer.MAX_VALUE --> 0.49999999976716936
+ nextInt == Integer.MIN_VALUE --> 0.5
+ nextInt == Integer.MAX_VALUE-1 --> 0.4999999995343387
+ nextInt == Integer.MIN_VALUE+1 --> 0.5000000002328306
+ nextInt == 1 --> 2.3283064365386963E-10
+ nextInt == -1 --> 0.9999999997671694
+ nextInt == 2 --> 4.6566128730773926E-10
+ nextInt == -2 --> 0.9999999995343387
+ */
+ }
+}
http://git-wip-us.apache.org/repos/asf/mahout/blob/545648f6/core/src/main/java/org/apache/mahout/math/jet/random/engine/package-info.java
----------------------------------------------------------------------
diff --git a/core/src/main/java/org/apache/mahout/math/jet/random/engine/package-info.java b/core/src/main/java/org/apache/mahout/math/jet/random/engine/package-info.java
new file mode 100644
index 0000000..e092010
--- /dev/null
+++ b/core/src/main/java/org/apache/mahout/math/jet/random/engine/package-info.java
@@ -0,0 +1,7 @@
+/**
+ * Engines generating strong uniformly distributed pseudo-random numbers;
+ * Needed by all JET probability distributions since they rely on uniform random numbers to generate random
+ * numbers from their own distribution.
+ * Thus, the classes of this package are at the core of computational statistics, simulations, Monte Carlo methods, etc.
+ */
+package org.apache.mahout.math.jet.random.engine;