Java Code Examples for org.apache.commons.math3.util.FastMath#random()

/**
 * Run with 77 random values, assuming that the outcome has a mean of 0 and a standard deviation of 1 with a
 * precision of 1E-10.
 */

@Test
public void testNormalize2() {
    // create an sample with 77 values    
    int length = 77;
    double sample[] = new double[length];
    for (int i = 0; i < length; i++) {
        sample[i] = FastMath.random();
    }
    // normalize this sample
    double standardizedSample[] = StatUtils.normalize(sample);

    DescriptiveStatistics stats = new DescriptiveStatistics();
    // Add the data from the array
    for (int i = 0; i < length; i++) {
        stats.addValue(standardizedSample[i]);
    }
    // the calculations do have a limited precision    
    double distance = 1E-10;
    // check the mean an standard deviation
    Assert.assertEquals(0.0, stats.getMean(), distance);
    Assert.assertEquals(1.0, stats.getStandardDeviation(), distance);

}
 

/**
 * Run with 77 random values, assuming that the outcome has a mean of 0 and a standard deviation of 1 with a
 * precision of 1E-10.
 */

@Test
public void testNormalize2() {
    // create an sample with 77 values    
    int length = 77;
    double sample[] = new double[length];
    for (int i = 0; i < length; i++) {
        sample[i] = FastMath.random();
    }
    // normalize this sample
    double standardizedSample[] = StatUtils.normalize(sample);

    DescriptiveStatistics stats = new DescriptiveStatistics();
    // Add the data from the array
    for (int i = 0; i < length; i++) {
        stats.addValue(standardizedSample[i]);
    }
    // the calculations do have a limited precision    
    double distance = 1E-10;
    // check the mean an standard deviation
    Assert.assertEquals(0.0, stats.getMean(), distance);
    Assert.assertEquals(1.0, stats.getStandardDeviation(), distance);

}
 

@Test
public void testSinc() {
    BivariateFunction div = new Divide();
    UnivariateFunction sin = new Sin();
    UnivariateFunction id = new Identity();
    UnivariateFunction sinc1 = FunctionUtils.combine(div, sin, id);
    UnivariateFunction sinc2 = new Sinc();

    for (int i = 0; i < 10; i++) {
        double x = FastMath.random();
        Assert.assertEquals(sinc1.value(x), sinc2.value(x), EPS);
    }
}
 

private void generateSineData(double[] xval, double[] yval, double xnoise, double ynoise) {
    double dx = 2 * FastMath.PI / xval.length;
    double x = 0;
    for(int i = 0; i < xval.length; ++i) {
        xval[i] = x;
        yval[i] = FastMath.sin(x) + (2 * FastMath.random() - 1) * ynoise;
        x += dx * (1 + (2 * FastMath.random() - 1) * xnoise);
    }
}
 

private void generateSineData(double[] xval, double[] yval, double xnoise, double ynoise) {
    double dx = 2 * FastMath.PI / xval.length;
    double x = 0;
    for(int i = 0; i < xval.length; ++i) {
        xval[i] = x;
        yval[i] = FastMath.sin(x) + (2 * FastMath.random() - 1) * ynoise;
        x += dx * (1 + (2 * FastMath.random() - 1) * xnoise);
    }
}
 

private void generateSineData(double[] xval, double[] yval, double xnoise, double ynoise) {
    double dx = 2 * FastMath.PI / xval.length;
    double x = 0;
    for(int i = 0; i < xval.length; ++i) {
        xval[i] = x;
        yval[i] = FastMath.sin(x) + (2 * FastMath.random() - 1) * ynoise;
        x += dx * (1 + (2 * FastMath.random() - 1) * xnoise);
    }
}
 

@Test
public void testFixingArguments() {
    UnivariateFunction scaler = FunctionUtils.fix1stArgument(new Multiply(), 10);
    Assert.assertEquals(1.23456, scaler.value(0.123456), EPS);

    UnivariateFunction pow1 = new Power(2);
    UnivariateFunction pow2 = FunctionUtils.fix2ndArgument(new Pow(), 2);

    for (int i = 0; i < 10; i++) {
        double x = FastMath.random() * 10;
        Assert.assertEquals(pow1.value(x), pow2.value(x), 0);
    }
}
 

private void generateSineData(double[] xval, double[] yval, double xnoise, double ynoise) {
    double dx = 2 * FastMath.PI / xval.length;
    double x = 0;
    for(int i = 0; i < xval.length; ++i) {
        xval[i] = x;
        yval[i] = FastMath.sin(x) + (2 * FastMath.random() - 1) * ynoise;
        x += dx * (1 + (2 * FastMath.random() - 1) * xnoise);
    }
}
 

private void generateSineData(double[] xval, double[] yval, double xnoise, double ynoise) {
    double dx = 2 * FastMath.PI / xval.length;
    double x = 0;
    for(int i = 0; i < xval.length; ++i) {
        xval[i] = x;
        yval[i] = FastMath.sin(x) + (2 * FastMath.random() - 1) * ynoise;
        x += dx * (1 + (2 * FastMath.random() - 1) * xnoise);
    }
}
 

private void generateSineData(double[] xval, double[] yval, double xnoise, double ynoise) {
    double dx = 2 * FastMath.PI / xval.length;
    double x = 0;
    for(int i = 0; i < xval.length; ++i) {
        xval[i] = x;
        yval[i] = FastMath.sin(x) + (2 * FastMath.random() - 1) * ynoise;
        x += dx * (1 + (2 * FastMath.random() - 1) * xnoise);
    }
}
 

private void generateSineData(double[] xval, double[] yval, double xnoise, double ynoise) {
    double dx = 2 * FastMath.PI / xval.length;
    double x = 0;
    for(int i = 0; i < xval.length; ++i) {
        xval[i] = x;
        yval[i] = FastMath.sin(x) + (2 * FastMath.random() - 1) * ynoise;
        x += dx * (1 + (2 * FastMath.random() - 1) * xnoise);
    }
}
 

@Test
public void testFixingArguments() {
    UnivariateFunction scaler = FunctionUtils.fix1stArgument(new Multiply(), 10);
    Assert.assertEquals(1.23456, scaler.value(0.123456), EPS);

    UnivariateFunction pow1 = new Power(2);
    UnivariateFunction pow2 = FunctionUtils.fix2ndArgument(new Pow(), 2);

    for (int i = 0; i < 10; i++) {
        double x = FastMath.random() * 10;
        Assert.assertEquals(pow1.value(x), pow2.value(x), 0);
    }
}
 

private void generateSineData(double[] xval, double[] yval, double xnoise, double ynoise) {
    double dx = 2 * FastMath.PI / xval.length;
    double x = 0;
    for(int i = 0; i < xval.length; ++i) {
        xval[i] = x;
        yval[i] = FastMath.sin(x) + (2 * FastMath.random() - 1) * ynoise;
        x += dx * (1 + (2 * FastMath.random() - 1) * xnoise);
    }
}
 

/**
 * {@inheritDoc}
 */
@Override
protected double doSolve()
    throws TooManyEvaluationsException,
           NumberIsTooLargeException,
           NoBracketingException {
    final double min = getMin();
    final double max = getMax();

    verifyInterval(min, max);

    final double relativeAccuracy = getRelativeAccuracy();
    final double absoluteAccuracy = getAbsoluteAccuracy();
    final double functionValueAccuracy = getFunctionValueAccuracy();

    // x2 is the last root approximation
    // x is the new approximation and new x2 for next round
    // x0 < x1 < x2 does not hold here

    double x0 = min;
    double y0 = computeObjectiveValue(x0);
    if (FastMath.abs(y0) < functionValueAccuracy) {
        return x0;
    }
    double x1 = max;
    double y1 = computeObjectiveValue(x1);
    if (FastMath.abs(y1) < functionValueAccuracy) {
        return x1;
    }

    if(y0 * y1 > 0) {
        throw new NoBracketingException(x0, x1, y0, y1);
    }

    double x2 = 0.5 * (x0 + x1);
    double y2 = computeObjectiveValue(x2);

    double oldx = Double.POSITIVE_INFINITY;
    while (true) {
        // quadratic interpolation through x0, x1, x2
        final double q = (x2 - x1) / (x1 - x0);
        final double a = q * (y2 - (1 + q) * y1 + q * y0);
        final double b = (2 * q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0;
        final double c = (1 + q) * y2;
        final double delta = b * b - 4 * a * c;
        double x;
        final double denominator;
        if (delta >= 0.0) {
            // choose a denominator larger in magnitude
            double dplus = b + FastMath.sqrt(delta);
            double dminus = b - FastMath.sqrt(delta);
            denominator = FastMath.abs(dplus) > FastMath.abs(dminus) ? dplus : dminus;
        } else {
            // take the modulus of (B +/- FastMath.sqrt(delta))
            denominator = FastMath.sqrt(b * b - delta);
        }
        if (denominator != 0) {
            x = x2 - 2.0 * c * (x2 - x1) / denominator;
            // perturb x if it exactly coincides with x1 or x2
            // the equality tests here are intentional
            while (x == x1 || x == x2) {
                x += absoluteAccuracy;
            }
        } else {
            // extremely rare case, get a random number to skip it
            x = min + FastMath.random() * (max - min);
            oldx = Double.POSITIVE_INFINITY;
        }
        final double y = computeObjectiveValue(x);

        // check for convergence
        final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy);
        if (FastMath.abs(x - oldx) <= tolerance ||
            FastMath.abs(y) <= functionValueAccuracy) {
            return x;
        }

        // prepare the next iteration
        x0 = x1;
        y0 = y1;
        x1 = x2;
        y1 = y2;
        x2 = x;
        y2 = y;
        oldx = x;
    }
}
 

/**
 * {@inheritDoc}
 */
@Override
protected double doSolve() {
    final double min = getMin();
    final double max = getMax();

    verifyInterval(min, max);

    final double relativeAccuracy = getRelativeAccuracy();
    final double absoluteAccuracy = getAbsoluteAccuracy();
    final double functionValueAccuracy = getFunctionValueAccuracy();

    // x2 is the last root approximation
    // x is the new approximation and new x2 for next round
    // x0 < x1 < x2 does not hold here

    double x0 = min;
    double y0 = computeObjectiveValue(x0);
    if (FastMath.abs(y0) < functionValueAccuracy) {
        return x0;
    }
    double x1 = max;
    double y1 = computeObjectiveValue(x1);
    if (FastMath.abs(y1) < functionValueAccuracy) {
        return x1;
    }

    if(y0 * y1 > 0) {
        throw new NoBracketingException(x0, x1, y0, y1);
    }

    double x2 = 0.5 * (x0 + x1);
    double y2 = computeObjectiveValue(x2);

    double oldx = Double.POSITIVE_INFINITY;
    while (true) {
        // quadratic interpolation through x0, x1, x2
        final double q = (x2 - x1) / (x1 - x0);
        final double a = q * (y2 - (1 + q) * y1 + q * y0);
        final double b = (2 * q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0;
        final double c = (1 + q) * y2;
        final double delta = b * b - 4 * a * c;
        double x;
        final double denominator;
        if (delta >= 0.0) {
            // choose a denominator larger in magnitude
            double dplus = b + FastMath.sqrt(delta);
            double dminus = b - FastMath.sqrt(delta);
            denominator = FastMath.abs(dplus) > FastMath.abs(dminus) ? dplus : dminus;
        } else {
            // take the modulus of (B +/- FastMath.sqrt(delta))
            denominator = FastMath.sqrt(b * b - delta);
        }
        if (denominator != 0) {
            x = x2 - 2.0 * c * (x2 - x1) / denominator;
            // perturb x if it exactly coincides with x1 or x2
            // the equality tests here are intentional
            while (x == x1 || x == x2) {
                x += absoluteAccuracy;
            }
        } else {
            // extremely rare case, get a random number to skip it
            x = min + FastMath.random() * (max - min);
            oldx = Double.POSITIVE_INFINITY;
        }
        final double y = computeObjectiveValue(x);

        // check for convergence
        final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy);
        if (FastMath.abs(x - oldx) <= tolerance ||
            FastMath.abs(y) <= functionValueAccuracy) {
            return x;
        }

        // prepare the next iteration
        x0 = x1;
        y0 = y1;
        x1 = x2;
        y1 = y2;
        x2 = x;
        y2 = y;
        oldx = x;
    }
}
 

/**
 * {@inheritDoc}
 */
@Override
protected double doSolve()
    throws TooManyEvaluationsException,
           NumberIsTooLargeException,
           NoBracketingException {
    final double min = getMin();
    final double max = getMax();

    verifyInterval(min, max);

    final double relativeAccuracy = getRelativeAccuracy();
    final double absoluteAccuracy = getAbsoluteAccuracy();
    final double functionValueAccuracy = getFunctionValueAccuracy();

    // x2 is the last root approximation
    // x is the new approximation and new x2 for next round
    // x0 < x1 < x2 does not hold here

    double x0 = min;
    double y0 = computeObjectiveValue(x0);
    if (FastMath.abs(y0) < functionValueAccuracy) {
        return x0;
    }
    double x1 = max;
    double y1 = computeObjectiveValue(x1);
    if (FastMath.abs(y1) < functionValueAccuracy) {
        return x1;
    }

    if(y0 * y1 > 0) {
        throw new NoBracketingException(x0, x1, y0, y1);
    }

    double x2 = 0.5 * (x0 + x1);
    double y2 = computeObjectiveValue(x2);

    double oldx = Double.POSITIVE_INFINITY;
    while (true) {
        // quadratic interpolation through x0, x1, x2
        final double q = (x2 - x1) / (x1 - x0);
        final double a = q * (y2 - (1 + q) * y1 + q * y0);
        final double b = (2 * q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0;
        final double c = (1 + q) * y2;
        final double delta = b * b - 4 * a * c;
        double x;
        final double denominator;
        if (delta >= 0.0) {
            // choose a denominator larger in magnitude
            double dplus = b + FastMath.sqrt(delta);
            double dminus = b - FastMath.sqrt(delta);
            denominator = FastMath.abs(dplus) > FastMath.abs(dminus) ? dplus : dminus;
        } else {
            // take the modulus of (B +/- FastMath.sqrt(delta))
            denominator = FastMath.sqrt(b * b - delta);
        }
        if (denominator != 0) {
            x = x2 - 2.0 * c * (x2 - x1) / denominator;
            // perturb x if it exactly coincides with x1 or x2
            // the equality tests here are intentional
            while (x == x1 || x == x2) {
                x += absoluteAccuracy;
            }
        } else {
            // extremely rare case, get a random number to skip it
            x = min + FastMath.random() * (max - min);
            oldx = Double.POSITIVE_INFINITY;
        }
        final double y = computeObjectiveValue(x);

        // check for convergence
        final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy);
        if (FastMath.abs(x - oldx) <= tolerance ||
            FastMath.abs(y) <= functionValueAccuracy) {
            return x;
        }

        // prepare the next iteration
        x0 = x1;
        y0 = y1;
        x1 = x2;
        y1 = y2;
        x2 = x;
        y2 = y;
        oldx = x;
    }
}
 

/**
 * {@inheritDoc}
 */
@Override
protected double doSolve()
    throws TooManyEvaluationsException,
           NumberIsTooLargeException,
           NoBracketingException {
    final double min = getMin();
    final double max = getMax();

    verifyInterval(min, max);

    final double relativeAccuracy = getRelativeAccuracy();
    final double absoluteAccuracy = getAbsoluteAccuracy();
    final double functionValueAccuracy = getFunctionValueAccuracy();

    // x2 is the last root approximation
    // x is the new approximation and new x2 for next round
    // x0 < x1 < x2 does not hold here

    double x0 = min;
    double y0 = computeObjectiveValue(x0);
    if (FastMath.abs(y0) < functionValueAccuracy) {
        return x0;
    }
    double x1 = max;
    double y1 = computeObjectiveValue(x1);
    if (FastMath.abs(y1) < functionValueAccuracy) {
        return x1;
    }

    if(y0 * y1 > 0) {
        throw new NoBracketingException(x0, x1, y0, y1);
    }

    double x2 = 0.5 * (x0 + x1);
    double y2 = computeObjectiveValue(x2);

    double oldx = Double.POSITIVE_INFINITY;
    while (true) {
        // quadratic interpolation through x0, x1, x2
        final double q = (x2 - x1) / (x1 - x0);
        final double a = q * (y2 - (1 + q) * y1 + q * y0);
        final double b = (2 * q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0;
        final double c = (1 + q) * y2;
        final double delta = b * b - 4 * a * c;
        double x;
        final double denominator;
        if (delta >= 0.0) {
            // choose a denominator larger in magnitude
            double dplus = b + FastMath.sqrt(delta);
            double dminus = b - FastMath.sqrt(delta);
            denominator = FastMath.abs(dplus) > FastMath.abs(dminus) ? dplus : dminus;
        } else {
            // take the modulus of (B +/- FastMath.sqrt(delta))
            denominator = FastMath.sqrt(b * b - delta);
        }
        if (denominator != 0) {
            x = x2 - 2.0 * c * (x2 - x1) / denominator;
            // perturb x if it exactly coincides with x1 or x2
            // the equality tests here are intentional
            while (x == x1 || x == x2) {
                x += absoluteAccuracy;
            }
        } else {
            // extremely rare case, get a random number to skip it
            x = min + FastMath.random() * (max - min);
            oldx = Double.POSITIVE_INFINITY;
        }
        final double y = computeObjectiveValue(x);

        // check for convergence
        final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy);
        if (FastMath.abs(x - oldx) <= tolerance ||
            FastMath.abs(y) <= functionValueAccuracy) {
            return x;
        }

        // prepare the next iteration
        x0 = x1;
        y0 = y1;
        x1 = x2;
        y1 = y2;
        x2 = x;
        y2 = y;
        oldx = x;
    }
}
 

/**
 * {@inheritDoc}
 */
@Override
protected double doSolve()
    throws TooManyEvaluationsException,
           NumberIsTooLargeException,
           NoBracketingException {
    final double min = getMin();
    final double max = getMax();

    verifyInterval(min, max);

    final double relativeAccuracy = getRelativeAccuracy();
    final double absoluteAccuracy = getAbsoluteAccuracy();
    final double functionValueAccuracy = getFunctionValueAccuracy();

    // x2 is the last root approximation
    // x is the new approximation and new x2 for next round
    // x0 < x1 < x2 does not hold here

    double x0 = min;
    double y0 = computeObjectiveValue(x0);
    if (FastMath.abs(y0) < functionValueAccuracy) {
        return x0;
    }
    double x1 = max;
    double y1 = computeObjectiveValue(x1);
    if (FastMath.abs(y1) < functionValueAccuracy) {
        return x1;
    }

    if(y0 * y1 > 0) {
        throw new NoBracketingException(x0, x1, y0, y1);
    }

    double x2 = 0.5 * (x0 + x1);
    double y2 = computeObjectiveValue(x2);

    double oldx = Double.POSITIVE_INFINITY;
    while (true) {
        // quadratic interpolation through x0, x1, x2
        final double q = (x2 - x1) / (x1 - x0);
        final double a = q * (y2 - (1 + q) * y1 + q * y0);
        final double b = (2 * q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0;
        final double c = (1 + q) * y2;
        final double delta = b * b - 4 * a * c;
        double x;
        final double denominator;
        if (delta >= 0.0) {
            // choose a denominator larger in magnitude
            double dplus = b + FastMath.sqrt(delta);
            double dminus = b - FastMath.sqrt(delta);
            denominator = FastMath.abs(dplus) > FastMath.abs(dminus) ? dplus : dminus;
        } else {
            // take the modulus of (B +/- FastMath.sqrt(delta))
            denominator = FastMath.sqrt(b * b - delta);
        }
        if (denominator != 0) {
            x = x2 - 2.0 * c * (x2 - x1) / denominator;
            // perturb x if it exactly coincides with x1 or x2
            // the equality tests here are intentional
            while (x == x1 || x == x2) {
                x += absoluteAccuracy;
            }
        } else {
            // extremely rare case, get a random number to skip it
            x = min + FastMath.random() * (max - min);
            oldx = Double.POSITIVE_INFINITY;
        }
        final double y = computeObjectiveValue(x);

        // check for convergence
        final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy);
        if (FastMath.abs(x - oldx) <= tolerance ||
            FastMath.abs(y) <= functionValueAccuracy) {
            return x;
        }

        // prepare the next iteration
        x0 = x1;
        y0 = y1;
        x1 = x2;
        y1 = y2;
        x2 = x;
        y2 = y;
        oldx = x;
    }
}
 

/**
 * {@inheritDoc}
 */
@Override
protected double doSolve()
    throws TooManyEvaluationsException,
           NumberIsTooLargeException,
           NoBracketingException {
    final double min = getMin();
    final double max = getMax();

    verifyInterval(min, max);

    final double relativeAccuracy = getRelativeAccuracy();
    final double absoluteAccuracy = getAbsoluteAccuracy();
    final double functionValueAccuracy = getFunctionValueAccuracy();

    // x2 is the last root approximation
    // x is the new approximation and new x2 for next round
    // x0 < x1 < x2 does not hold here

    double x0 = min;
    double y0 = computeObjectiveValue(x0);
    if (FastMath.abs(y0) < functionValueAccuracy) {
        return x0;
    }
    double x1 = max;
    double y1 = computeObjectiveValue(x1);
    if (FastMath.abs(y1) < functionValueAccuracy) {
        return x1;
    }

    if(y0 * y1 > 0) {
        throw new NoBracketingException(x0, x1, y0, y1);
    }

    double x2 = 0.5 * (x0 + x1);
    double y2 = computeObjectiveValue(x2);

    double oldx = Double.POSITIVE_INFINITY;
    while (true) {
        // quadratic interpolation through x0, x1, x2
        final double q = (x2 - x1) / (x1 - x0);
        final double a = q * (y2 - (1 + q) * y1 + q * y0);
        final double b = (2 * q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0;
        final double c = (1 + q) * y2;
        final double delta = b * b - 4 * a * c;
        double x;
        final double denominator;
        if (delta >= 0.0) {
            // choose a denominator larger in magnitude
            double dplus = b + FastMath.sqrt(delta);
            double dminus = b - FastMath.sqrt(delta);
            denominator = FastMath.abs(dplus) > FastMath.abs(dminus) ? dplus : dminus;
        } else {
            // take the modulus of (B +/- FastMath.sqrt(delta))
            denominator = FastMath.sqrt(b * b - delta);
        }
        if (denominator != 0) {
            x = x2 - 2.0 * c * (x2 - x1) / denominator;
            // perturb x if it exactly coincides with x1 or x2
            // the equality tests here are intentional
            while (x == x1 || x == x2) {
                x += absoluteAccuracy;
            }
        } else {
            // extremely rare case, get a random number to skip it
            x = min + FastMath.random() * (max - min);
            oldx = Double.POSITIVE_INFINITY;
        }
        final double y = computeObjectiveValue(x);

        // check for convergence
        final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy);
        if (FastMath.abs(x - oldx) <= tolerance ||
            FastMath.abs(y) <= functionValueAccuracy) {
            return x;
        }

        // prepare the next iteration
        x0 = x1;
        y0 = y1;
        x1 = x2;
        y1 = y2;
        x2 = x;
        y2 = y;
        oldx = x;
    }
}
 

/**
 * {@inheritDoc}
 */
@Override
protected double doSolve()
    throws TooManyEvaluationsException,
           NumberIsTooLargeException,
           NoBracketingException {
    final double min = getMin();
    final double max = getMax();

    verifyInterval(min, max);

    final double relativeAccuracy = getRelativeAccuracy();
    final double absoluteAccuracy = getAbsoluteAccuracy();
    final double functionValueAccuracy = getFunctionValueAccuracy();

    // x2 is the last root approximation
    // x is the new approximation and new x2 for next round
    // x0 < x1 < x2 does not hold here

    double x0 = min;
    double y0 = computeObjectiveValue(x0);
    if (FastMath.abs(y0) < functionValueAccuracy) {
        return x0;
    }
    double x1 = max;
    double y1 = computeObjectiveValue(x1);
    if (FastMath.abs(y1) < functionValueAccuracy) {
        return x1;
    }

    if(y0 * y1 > 0) {
        throw new NoBracketingException(x0, x1, y0, y1);
    }

    double x2 = 0.5 * (x0 + x1);
    double y2 = computeObjectiveValue(x2);

    double oldx = Double.POSITIVE_INFINITY;
    while (true) {
        // quadratic interpolation through x0, x1, x2
        final double q = (x2 - x1) / (x1 - x0);
        final double a = q * (y2 - (1 + q) * y1 + q * y0);
        final double b = (2 * q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0;
        final double c = (1 + q) * y2;
        final double delta = b * b - 4 * a * c;
        double x;
        final double denominator;
        if (delta >= 0.0) {
            // choose a denominator larger in magnitude
            double dplus = b + FastMath.sqrt(delta);
            double dminus = b - FastMath.sqrt(delta);
            denominator = FastMath.abs(dplus) > FastMath.abs(dminus) ? dplus : dminus;
        } else {
            // take the modulus of (B +/- FastMath.sqrt(delta))
            denominator = FastMath.sqrt(b * b - delta);
        }
        if (denominator != 0) {
            x = x2 - 2.0 * c * (x2 - x1) / denominator;
            // perturb x if it exactly coincides with x1 or x2
            // the equality tests here are intentional
            while (x == x1 || x == x2) {
                x += absoluteAccuracy;
            }
        } else {
            // extremely rare case, get a random number to skip it
            x = min + FastMath.random() * (max - min);
            oldx = Double.POSITIVE_INFINITY;
        }
        final double y = computeObjectiveValue(x);

        // check for convergence
        final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy);
        if (FastMath.abs(x - oldx) <= tolerance ||
            FastMath.abs(y) <= functionValueAccuracy) {
            return x;
        }

        // prepare the next iteration
        x0 = x1;
        y0 = y1;
        x1 = x2;
        y1 = y2;
        x2 = x;
        y2 = y;
        oldx = x;
    }
}