Java Code Examples for org.apache.commons.math3.exception.util.LocalizedFormats#ZERO_DENOMINATOR

/** Finds r such that a = q b + r with 0 <= r < b if b > 0 and b < r <= 0 if b > 0.
 * <p>
 * This methods returns the same value as integer modulo when
 * a and b are same signs, but returns a different value when
 * they are opposite (i.e. q is negative).
 * </p>
 * @param a dividend
 * @param b divisor
 * @return r such that a = q b + r with 0 <= r < b if b > 0 and b < r <= 0 if b > 0
 * @exception MathArithmeticException if b == 0
 * @see #floorDiv(int, int)
 * @since 3.4
 */
public static int floorMod(final int a, final int b) throws MathArithmeticException {

    if (b == 0) {
        throw new MathArithmeticException(LocalizedFormats.ZERO_DENOMINATOR);
    }

    final int m = a % b;
    if ((a ^ b) >= 0 || m == 0) {
        // a an b have same sign, or division is exact
        return m;
    } else {
        // a and b have opposite signs and division is not exact
        return b + m;
    }

}
 

/**Solve  a  system of composed of a Lower Triangular Matrix
 * {@link RealMatrix}.
 * <p>
 * This method is called to solve systems of equations which are
 * of the lower triangular form. The matrix {@link RealMatrix}
 * is assumed, though not checked, to be in lower triangular form.
 * The vector {@link RealVector} is overwritten with the solution.
 * The matrix is checked that it is square and its dimensions match
 * the length of the vector.
 * </p>
 * @param rm RealMatrix which is lower triangular
 * @param b  RealVector this is overwritten
 * @throws DimensionMismatchException if the matrix and vector are not
 * conformable
 * @throws NonSquareMatrixException if the matrix {@code rm} is not square
 * @throws MathArithmeticException if the absolute value of one of the diagonal
 * coefficient of {@code rm} is lower than {@link Precision#SAFE_MIN}
 */
public static void solveLowerTriangularSystem(RealMatrix rm, RealVector b)
    throws DimensionMismatchException, MathArithmeticException,
    NonSquareMatrixException {
    if ((rm == null) || (b == null) || ( rm.getRowDimension() != b.getDimension())) {
        throw new DimensionMismatchException(
                (rm == null) ? 0 : rm.getRowDimension(),
                (b == null) ? 0 : b.getDimension());
    }
    if( rm.getColumnDimension() != rm.getRowDimension() ){
        throw new NonSquareMatrixException(rm.getRowDimension(),
                                           rm.getColumnDimension());
    }
    int rows = rm.getRowDimension();
    for( int i = 0 ; i < rows ; i++ ){
        double diag = rm.getEntry(i, i);
        if( FastMath.abs(diag) < Precision.SAFE_MIN ){
            throw new MathArithmeticException(LocalizedFormats.ZERO_DENOMINATOR);
        }
        double bi = b.getEntry(i)/diag;
        b.setEntry(i,  bi );
        for( int j = i+1; j< rows; j++ ){
            b.setEntry(j, b.getEntry(j)-bi*rm.getEntry(j,i)  );
        }
    }
}
 

/** Finds q such that a = q b + r with 0 <= r < b if b > 0 and b < r <= 0 if b > 0.
 * <p>
 * This methods returns the same value as integer division when
 * a and b are same signs, but returns a different value when
 * they are opposite (i.e. q is negative).
 * </p>
 * @param a dividend
 * @param b divisor
 * @return q such that a = q b + r with 0 <= r < b if b > 0 and b < r <= 0 if b > 0
 * @exception MathArithmeticException if b == 0
 * @see #floorMod(long, long)
 * @since 3.4
 */
public static long floorDiv(final long a, final long b) throws MathArithmeticException {

    if (b == 0l) {
        throw new MathArithmeticException(LocalizedFormats.ZERO_DENOMINATOR);
    }

    final long m = a % b;
    if ((a ^ b) >= 0l || m == 0l) {
        // a an b have same sign, or division is exact
        return a / b;
    } else {
        // a and b have opposite signs and division is not exact
        return (a / b) - 1l;
    }

}
 

/**Solve  a  system of composed of a Lower Triangular Matrix
 * {@link RealMatrix}.
 * <p>
 * This method is called to solve systems of equations which are
 * of the lower triangular form. The matrix {@link RealMatrix}
 * is assumed, though not checked, to be in lower triangular form.
 * The vector {@link RealVector} is overwritten with the solution.
 * The matrix is checked that it is square and its dimensions match
 * the length of the vector.
 * </p>
 * @param rm RealMatrix which is lower triangular
 * @param b  RealVector this is overwritten
 * @throws DimensionMismatchException if the matrix and vector are not
 * conformable
 * @throws NonSquareMatrixException if the matrix {@code rm} is not square
 * @throws MathArithmeticException if the absolute value of one of the diagonal
 * coefficient of {@code rm} is lower than {@link Precision#SAFE_MIN}
 */
public static void solveLowerTriangularSystem(RealMatrix rm, RealVector b)
    throws DimensionMismatchException, MathArithmeticException,
    NonSquareMatrixException {
    if ((rm == null) || (b == null) || ( rm.getRowDimension() != b.getDimension())) {
        throw new DimensionMismatchException(
                (rm == null) ? 0 : rm.getRowDimension(),
                (b == null) ? 0 : b.getDimension());
    }
    if( rm.getColumnDimension() != rm.getRowDimension() ){
        throw new NonSquareMatrixException(rm.getRowDimension(),
                                           rm.getColumnDimension());
    }
    int rows = rm.getRowDimension();
    for( int i = 0 ; i < rows ; i++ ){
        double diag = rm.getEntry(i, i);
        if( FastMath.abs(diag) < Precision.SAFE_MIN ){
            throw new MathArithmeticException(LocalizedFormats.ZERO_DENOMINATOR);
        }
        double bi = b.getEntry(i)/diag;
        b.setEntry(i,  bi );
        for( int j = i+1; j< rows; j++ ){
            b.setEntry(j, b.getEntry(j)-bi*rm.getEntry(j,i)  );
        }
    }
}
 

/** Solver a  system composed  of an Upper Triangular Matrix
 * {@link RealMatrix}.
 * <p>
 * This method is called to solve systems of equations which are
 * of the lower triangular form. The matrix {@link RealMatrix}
 * is assumed, though not checked, to be in upper triangular form.
 * The vector {@link RealVector} is overwritten with the solution.
 * The matrix is checked that it is square and its dimensions match
 * the length of the vector.
 * </p>
 * @param rm RealMatrix which is upper triangular
 * @param b  RealVector this is overwritten
 * @exception IllegalArgumentException if the matrix and vector are not conformable
 * @exception ArithmeticException there is a zero or near zero on the diagonal of rm
 */
public static void solveUpperTriangularSystem( RealMatrix rm, RealVector b){
    if ((rm == null) || (b == null) || ( rm.getRowDimension() != b.getDimension())) {
        throw new MathIllegalArgumentException(LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE,
                (rm == null) ? 0 : rm.getRowDimension(),
                (b == null) ? 0 : b.getDimension());
    }
    if( rm.getColumnDimension() != rm.getRowDimension() ){
        throw new MathIllegalArgumentException(LocalizedFormats.DIMENSIONS_MISMATCH_2x2,
                rm.getRowDimension(),rm.getRowDimension(),
                rm.getRowDimension(),rm.getColumnDimension());
    }
    int rows = rm.getRowDimension();
    for( int i = rows-1 ; i >-1 ; i-- ){
        double diag = rm.getEntry(i, i);
        if( FastMath.abs(diag) < Precision.SAFE_MIN ){
            throw new MathArithmeticException(LocalizedFormats.ZERO_DENOMINATOR);
        }
        double bi = b.getEntry(i)/diag;
        b.setEntry(i,  bi );
        for( int j = i-1; j>-1; j-- ){
            b.setEntry(j, b.getEntry(j)-bi*rm.getEntry(j,i)  );
        }
    }
}
 

/**
 * Create a {@link BigFraction} given the numerator and denominator as
 * {@code BigInteger}. The {@link BigFraction} is reduced to lowest terms.
 *
 * @param num the numerator, must not be {@code null}.
 * @param den the denominator, must not be {@code null}.
 * @throws ZeroException if the denominator is zero.
 * @throws NullArgumentException if either of the arguments is null
 */
public BigFraction(BigInteger num, BigInteger den) {
    MathUtils.checkNotNull(num, LocalizedFormats.NUMERATOR);
    MathUtils.checkNotNull(den, LocalizedFormats.DENOMINATOR);
    if (BigInteger.ZERO.equals(den)) {
        throw new ZeroException(LocalizedFormats.ZERO_DENOMINATOR);
    }
    if (BigInteger.ZERO.equals(num)) {
        numerator   = BigInteger.ZERO;
        denominator = BigInteger.ONE;
    } else {

        // reduce numerator and denominator by greatest common denominator
        final BigInteger gcd = num.gcd(den);
        if (BigInteger.ONE.compareTo(gcd) < 0) {
            num = num.divide(gcd);
            den = den.divide(gcd);
        }

        // move sign to numerator
        if (BigInteger.ZERO.compareTo(den) > 0) {
            num = num.negate();
            den = den.negate();
        }

        // store the values in the final fields
        numerator   = num;
        denominator = den;

    }
}
 

/**
 * Create a {@link BigFraction} given the numerator and denominator as
 * {@code BigInteger}. The {@link BigFraction} is reduced to lowest terms.
 *
 * @param num the numerator, must not be {@code null}.
 * @param den the denominator, must not be {@code null}.
 * @throws ZeroException if the denominator is zero.
 * @throws NullArgumentException if either of the arguments is null
 */
public BigFraction(BigInteger num, BigInteger den) {
    MathUtils.checkNotNull(num, LocalizedFormats.NUMERATOR);
    MathUtils.checkNotNull(den, LocalizedFormats.DENOMINATOR);
    if (BigInteger.ZERO.equals(den)) {
        throw new ZeroException(LocalizedFormats.ZERO_DENOMINATOR);
    }
    if (BigInteger.ZERO.equals(num)) {
        numerator   = BigInteger.ZERO;
        denominator = BigInteger.ONE;
    } else {

        // reduce numerator and denominator by greatest common denominator
        final BigInteger gcd = num.gcd(den);
        if (BigInteger.ONE.compareTo(gcd) < 0) {
            num = num.divide(gcd);
            den = den.divide(gcd);
        }

        // move sign to numerator
        if (BigInteger.ZERO.compareTo(den) > 0) {
            num = num.negate();
            den = den.negate();
        }

        // store the values in the final fields
        numerator   = num;
        denominator = den;

    }
}
 

/**
 * Create a {@link BigFraction} given the numerator and denominator as
 * {@code BigInteger}. The {@link BigFraction} is reduced to lowest terms.
 *
 * @param num the numerator, must not be {@code null}.
 * @param den the denominator, must not be {@code null}.
 * @throws ZeroException if the denominator is zero.
 * @throws NullArgumentException if either of the arguments is null
 */
public BigFraction(BigInteger num, BigInteger den) {
    MathUtils.checkNotNull(num, LocalizedFormats.NUMERATOR);
    MathUtils.checkNotNull(den, LocalizedFormats.DENOMINATOR);
    if (BigInteger.ZERO.equals(den)) {
        throw new ZeroException(LocalizedFormats.ZERO_DENOMINATOR);
    }
    if (BigInteger.ZERO.equals(num)) {
        numerator   = BigInteger.ZERO;
        denominator = BigInteger.ONE;
    } else {

        // reduce numerator and denominator by greatest common denominator
        final BigInteger gcd = num.gcd(den);
        if (BigInteger.ONE.compareTo(gcd) < 0) {
            num = num.divide(gcd);
            den = den.divide(gcd);
        }

        // move sign to numerator
        if (BigInteger.ZERO.compareTo(den) > 0) {
            num = num.negate();
            den = den.negate();
        }

        // store the values in the final fields
        numerator   = num;
        denominator = den;

    }
}
 

/**
 * Create a {@link BigFraction} given the numerator and denominator as
 * {@code BigInteger}. The {@link BigFraction} is reduced to lowest terms.
 *
 * @param num the numerator, must not be {@code null}.
 * @param den the denominator, must not be {@code null}.
 * @throws ZeroException if the denominator is zero.
 * @throws NullArgumentException if either of the arguments is null
 */
public BigFraction(BigInteger num, BigInteger den) {
    MathUtils.checkNotNull(num, LocalizedFormats.NUMERATOR);
    MathUtils.checkNotNull(den, LocalizedFormats.DENOMINATOR);
    if (BigInteger.ZERO.equals(den)) {
        throw new ZeroException(LocalizedFormats.ZERO_DENOMINATOR);
    }
    if (BigInteger.ZERO.equals(num)) {
        numerator   = BigInteger.ZERO;
        denominator = BigInteger.ONE;
    } else {

        // reduce numerator and denominator by greatest common denominator
        final BigInteger gcd = num.gcd(den);
        if (BigInteger.ONE.compareTo(gcd) < 0) {
            num = num.divide(gcd);
            den = den.divide(gcd);
        }

        // move sign to numerator
        if (BigInteger.ZERO.compareTo(den) > 0) {
            num = num.negate();
            den = den.negate();
        }

        // store the values in the final fields
        numerator   = num;
        denominator = den;

    }
}
 

/**
 * Estimate a first guess of the amplitude and angular frequency.
 * This method assumes that the {@link #sortObservations(WeightedObservedPoint[])} method
 * has been called previously.
 *
 * @param observations Observations, sorted w.r.t. abscissa.
 * @throws ZeroException if the abscissa range is zero.
 * @throws MathIllegalStateException when the guessing procedure cannot
 * produce sensible results.
 * @return the guessed amplitude (at index 0) and circular frequency
 * (at index 1).
 */
private double[] guessAOmega(WeightedObservedPoint[] observations) {
    final double[] aOmega = new double[2];

    // initialize the sums for the linear model between the two integrals
    double sx2 = 0;
    double sy2 = 0;
    double sxy = 0;
    double sxz = 0;
    double syz = 0;

    double currentX = observations[0].getX();
    double currentY = observations[0].getY();
    double f2Integral = 0;
    double fPrime2Integral = 0;
    final double startX = currentX;
    for (int i = 1; i < observations.length; ++i) {
        // one step forward
        final double previousX = currentX;
        final double previousY = currentY;
        currentX = observations[i].getX();
        currentY = observations[i].getY();

        // update the integrals of f<sup>2</sup> and f'<sup>2</sup>
        // considering a linear model for f (and therefore constant f')
        final double dx = currentX - previousX;
        final double dy = currentY - previousY;
        final double f2StepIntegral =
            dx * (previousY * previousY + previousY * currentY + currentY * currentY) / 3;
        final double fPrime2StepIntegral = dy * dy / dx;

        final double x = currentX - startX;
        f2Integral += f2StepIntegral;
        fPrime2Integral += fPrime2StepIntegral;

        sx2 += x * x;
        sy2 += f2Integral * f2Integral;
        sxy += x * f2Integral;
        sxz += x * fPrime2Integral;
        syz += f2Integral * fPrime2Integral;
    }

    // compute the amplitude and pulsation coefficients
    double c1 = sy2 * sxz - sxy * syz;
    double c2 = sxy * sxz - sx2 * syz;
    double c3 = sx2 * sy2 - sxy * sxy;
    if ((c1 / c2 < 0) || (c2 / c3 < 0)) {
        final int last = observations.length - 1;
        // Range of the observations, assuming that the
        // observations are sorted.
        final double xRange = observations[last].getX() - observations[0].getX();
        if (xRange == 0) {
            throw new ZeroException();
        }
        aOmega[1] = 2 * Math.PI / xRange;

        double yMin = Double.POSITIVE_INFINITY;
        double yMax = Double.NEGATIVE_INFINITY;
        for (int i = 1; i < observations.length; ++i) {
            final double y = observations[i].getY();
            if (y < yMin) {
                yMin = y;
            }
            if (y > yMax) {
                yMax = y;
            }
        }
        aOmega[0] = 0.5 * (yMax - yMin);
    } else {
        if (c2 == 0) {
            // In some ill-conditioned cases (cf. MATH-844), the guesser
            // procedure cannot produce sensible results.
            throw new MathIllegalStateException(LocalizedFormats.ZERO_DENOMINATOR);
        }

        aOmega[0] = FastMath.sqrt(c1 / c2);
        aOmega[1] = FastMath.sqrt(c2 / c3);
    }

    return aOmega;
}
 

/**
 * Estimate a first guess of the amplitude and angular frequency.
 * This method assumes that the {@link #sortObservations()} method
 * has been called previously.
 *
 * @param observations Observations, sorted w.r.t. abscissa.
 * @throws ZeroException if the abscissa range is zero.
 * @throws MathIllegalStateException when the guessing procedure cannot
 * produce sensible results.
 * @return the guessed amplitude (at index 0) and circular frequency
 * (at index 1).
 */
private double[] guessAOmega(WeightedObservedPoint[] observations) {
    final double[] aOmega = new double[2];

    // initialize the sums for the linear model between the two integrals
    double sx2 = 0;
    double sy2 = 0;
    double sxy = 0;
    double sxz = 0;
    double syz = 0;

    double currentX = observations[0].getX();
    double currentY = observations[0].getY();
    double f2Integral = 0;
    double fPrime2Integral = 0;
    final double startX = currentX;
    for (int i = 1; i < observations.length; ++i) {
        // one step forward
        final double previousX = currentX;
        final double previousY = currentY;
        currentX = observations[i].getX();
        currentY = observations[i].getY();

        // update the integrals of f<sup>2</sup> and f'<sup>2</sup>
        // considering a linear model for f (and therefore constant f')
        final double dx = currentX - previousX;
        final double dy = currentY - previousY;
        final double f2StepIntegral =
            dx * (previousY * previousY + previousY * currentY + currentY * currentY) / 3;
        final double fPrime2StepIntegral = dy * dy / dx;

        final double x = currentX - startX;
        f2Integral += f2StepIntegral;
        fPrime2Integral += fPrime2StepIntegral;

        sx2 += x * x;
        sy2 += f2Integral * f2Integral;
        sxy += x * f2Integral;
        sxz += x * fPrime2Integral;
        syz += f2Integral * fPrime2Integral;
    }

    // compute the amplitude and pulsation coefficients
    double c1 = sy2 * sxz - sxy * syz;
    double c2 = sxy * sxz - sx2 * syz;
    double c3 = sx2 * sy2 - sxy * sxy;
    if ((c1 / c2 < 0) || (c2 / c3 < 0)) {
        final int last = observations.length - 1;
        // Range of the observations, assuming that the
        // observations are sorted.
        final double xRange = observations[last].getX() - observations[0].getX();
        if (xRange == 0) {
            throw new ZeroException();
        }
        aOmega[1] = 2 * Math.PI / xRange;

        double yMin = Double.POSITIVE_INFINITY;
        double yMax = Double.NEGATIVE_INFINITY;
        for (int i = 1; i < observations.length; ++i) {
            final double y = observations[i].getY();
            if (y < yMin) {
                yMin = y;
            }
            if (y > yMax) {
                yMax = y;
            }
        }
        aOmega[0] = 0.5 * (yMax - yMin);
    } else {
        if (c2 == 0) {
            // In some ill-conditioned cases (cf. MATH-844), the guesser
            // procedure cannot produce sensible results.
            throw new MathIllegalStateException(LocalizedFormats.ZERO_DENOMINATOR);
        }

        aOmega[0] = FastMath.sqrt(c1 / c2);
        aOmega[1] = FastMath.sqrt(c2 / c3);
    }

    return aOmega;
}
 

/**
 * Estimate a first guess of the amplitude and angular frequency.
 *
 * @param observations Observations, sorted w.r.t. abscissa.
 * @throws ZeroException if the abscissa range is zero.
 * @throws MathIllegalStateException when the guessing procedure cannot
 * produce sensible results.
 * @return the guessed amplitude (at index 0) and circular frequency
 * (at index 1).
 */
private double[] guessAOmega(WeightedObservedPoint[] observations) {
    final double[] aOmega = new double[2];

    // initialize the sums for the linear model between the two integrals
    double sx2 = 0;
    double sy2 = 0;
    double sxy = 0;
    double sxz = 0;
    double syz = 0;

    double currentX = observations[0].getX();
    double currentY = observations[0].getY();
    double f2Integral = 0;
    double fPrime2Integral = 0;
    final double startX = currentX;
    for (int i = 1; i < observations.length; ++i) {
        // one step forward
        final double previousX = currentX;
        final double previousY = currentY;
        currentX = observations[i].getX();
        currentY = observations[i].getY();

        // update the integrals of f<sup>2</sup> and f'<sup>2</sup>
        // considering a linear model for f (and therefore constant f')
        final double dx = currentX - previousX;
        final double dy = currentY - previousY;
        final double f2StepIntegral =
            dx * (previousY * previousY + previousY * currentY + currentY * currentY) / 3;
        final double fPrime2StepIntegral = dy * dy / dx;

        final double x = currentX - startX;
        f2Integral += f2StepIntegral;
        fPrime2Integral += fPrime2StepIntegral;

        sx2 += x * x;
        sy2 += f2Integral * f2Integral;
        sxy += x * f2Integral;
        sxz += x * fPrime2Integral;
        syz += f2Integral * fPrime2Integral;
    }

    // compute the amplitude and pulsation coefficients
    double c1 = sy2 * sxz - sxy * syz;
    double c2 = sxy * sxz - sx2 * syz;
    double c3 = sx2 * sy2 - sxy * sxy;
    if ((c1 / c2 < 0) || (c2 / c3 < 0)) {
        final int last = observations.length - 1;
        // Range of the observations, assuming that the
        // observations are sorted.
        final double xRange = observations[last].getX() - observations[0].getX();
        if (xRange == 0) {
            throw new ZeroException();
        }
        aOmega[1] = 2 * Math.PI / xRange;

        double yMin = Double.POSITIVE_INFINITY;
        double yMax = Double.NEGATIVE_INFINITY;
        for (int i = 1; i < observations.length; ++i) {
            final double y = observations[i].getY();
            if (y < yMin) {
                yMin = y;
            }
            if (y > yMax) {
                yMax = y;
            }
        }
        aOmega[0] = 0.5 * (yMax - yMin);
    } else {
        if (c2 == 0) {
            // In some ill-conditioned cases (cf. MATH-844), the guesser
            // procedure cannot produce sensible results.
            throw new MathIllegalStateException(LocalizedFormats.ZERO_DENOMINATOR);
        }

        aOmega[0] = FastMath.sqrt(c1 / c2);
        aOmega[1] = FastMath.sqrt(c2 / c3);
    }

    return aOmega;
}
 

/**
 * Estimate a first guess of the amplitude and angular frequency.
 * This method assumes that the {@link #sortObservations(WeightedObservedPoint[])} method
 * has been called previously.
 *
 * @param observations Observations, sorted w.r.t. abscissa.
 * @throws ZeroException if the abscissa range is zero.
 * @throws MathIllegalStateException when the guessing procedure cannot
 * produce sensible results.
 * @return the guessed amplitude (at index 0) and circular frequency
 * (at index 1).
 */
private double[] guessAOmega(WeightedObservedPoint[] observations) {
    final double[] aOmega = new double[2];

    // initialize the sums for the linear model between the two integrals
    double sx2 = 0;
    double sy2 = 0;
    double sxy = 0;
    double sxz = 0;
    double syz = 0;

    double currentX = observations[0].getX();
    double currentY = observations[0].getY();
    double f2Integral = 0;
    double fPrime2Integral = 0;
    final double startX = currentX;
    for (int i = 1; i < observations.length; ++i) {
        // one step forward
        final double previousX = currentX;
        final double previousY = currentY;
        currentX = observations[i].getX();
        currentY = observations[i].getY();

        // update the integrals of f<sup>2</sup> and f'<sup>2</sup>
        // considering a linear model for f (and therefore constant f')
        final double dx = currentX - previousX;
        final double dy = currentY - previousY;
        final double f2StepIntegral =
            dx * (previousY * previousY + previousY * currentY + currentY * currentY) / 3;
        final double fPrime2StepIntegral = dy * dy / dx;

        final double x = currentX - startX;
        f2Integral += f2StepIntegral;
        fPrime2Integral += fPrime2StepIntegral;

        sx2 += x * x;
        sy2 += f2Integral * f2Integral;
        sxy += x * f2Integral;
        sxz += x * fPrime2Integral;
        syz += f2Integral * fPrime2Integral;
    }

    // compute the amplitude and pulsation coefficients
    double c1 = sy2 * sxz - sxy * syz;
    double c2 = sxy * sxz - sx2 * syz;
    double c3 = sx2 * sy2 - sxy * sxy;
    if ((c1 / c2 < 0) || (c2 / c3 < 0)) {
        final int last = observations.length - 1;
        // Range of the observations, assuming that the
        // observations are sorted.
        final double xRange = observations[last].getX() - observations[0].getX();
        if (xRange == 0) {
            throw new ZeroException();
        }
        aOmega[1] = 2 * Math.PI / xRange;

        double yMin = Double.POSITIVE_INFINITY;
        double yMax = Double.NEGATIVE_INFINITY;
        for (int i = 1; i < observations.length; ++i) {
            final double y = observations[i].getY();
            if (y < yMin) {
                yMin = y;
            }
            if (y > yMax) {
                yMax = y;
            }
        }
        aOmega[0] = 0.5 * (yMax - yMin);
    } else {
        if (c2 == 0) {
            // In some ill-conditioned cases (cf. MATH-844), the guesser
            // procedure cannot produce sensible results.
            throw new MathIllegalStateException(LocalizedFormats.ZERO_DENOMINATOR);
        }

        aOmega[0] = FastMath.sqrt(c1 / c2);
        aOmega[1] = FastMath.sqrt(c2 / c3);
    }

    return aOmega;
}
 

/**
 * Estimate a first guess of the amplitude and angular frequency.
 *
 * @param observations Observations, sorted w.r.t. abscissa.
 * @throws ZeroException if the abscissa range is zero.
 * @throws MathIllegalStateException when the guessing procedure cannot
 * produce sensible results.
 * @return the guessed amplitude (at index 0) and circular frequency
 * (at index 1).
 */
private double[] guessAOmega(WeightedObservedPoint[] observations) {
    final double[] aOmega = new double[2];

    // initialize the sums for the linear model between the two integrals
    double sx2 = 0;
    double sy2 = 0;
    double sxy = 0;
    double sxz = 0;
    double syz = 0;

    double currentX = observations[0].getX();
    double currentY = observations[0].getY();
    double f2Integral = 0;
    double fPrime2Integral = 0;
    final double startX = currentX;
    for (int i = 1; i < observations.length; ++i) {
        // one step forward
        final double previousX = currentX;
        final double previousY = currentY;
        currentX = observations[i].getX();
        currentY = observations[i].getY();

        // update the integrals of f<sup>2</sup> and f'<sup>2</sup>
        // considering a linear model for f (and therefore constant f')
        final double dx = currentX - previousX;
        final double dy = currentY - previousY;
        final double f2StepIntegral =
            dx * (previousY * previousY + previousY * currentY + currentY * currentY) / 3;
        final double fPrime2StepIntegral = dy * dy / dx;

        final double x = currentX - startX;
        f2Integral += f2StepIntegral;
        fPrime2Integral += fPrime2StepIntegral;

        sx2 += x * x;
        sy2 += f2Integral * f2Integral;
        sxy += x * f2Integral;
        sxz += x * fPrime2Integral;
        syz += f2Integral * fPrime2Integral;
    }

    // compute the amplitude and pulsation coefficients
    double c1 = sy2 * sxz - sxy * syz;
    double c2 = sxy * sxz - sx2 * syz;
    double c3 = sx2 * sy2 - sxy * sxy;
    if ((c1 / c2 < 0) || (c2 / c3 < 0)) {
        final int last = observations.length - 1;
        // Range of the observations, assuming that the
        // observations are sorted.
        final double xRange = observations[last].getX() - observations[0].getX();
        if (xRange == 0) {
            throw new ZeroException();
        }
        aOmega[1] = 2 * Math.PI / xRange;

        double yMin = Double.POSITIVE_INFINITY;
        double yMax = Double.NEGATIVE_INFINITY;
        for (int i = 1; i < observations.length; ++i) {
            final double y = observations[i].getY();
            if (y < yMin) {
                yMin = y;
            }
            if (y > yMax) {
                yMax = y;
            }
        }
        aOmega[0] = 0.5 * (yMax - yMin);
    } else {
        if (c2 == 0) {
            // In some ill-conditioned cases (cf. MATH-844), the guesser
            // procedure cannot produce sensible results.
            throw new MathIllegalStateException(LocalizedFormats.ZERO_DENOMINATOR);
        }

        aOmega[0] = FastMath.sqrt(c1 / c2);
        aOmega[1] = FastMath.sqrt(c2 / c3);
    }

    return aOmega;
}
 

/**
 * Estimate a first guess of the amplitude and angular frequency.
 * This method assumes that the {@link #sortObservations()} method
 * has been called previously.
 *
 * @param observations Observations, sorted w.r.t. abscissa.
 * @throws ZeroException if the abscissa range is zero.
 * @throws MathIllegalStateException when the guessing procedure cannot
 * produce sensible results.
 * @return the guessed amplitude (at index 0) and circular frequency
 * (at index 1).
 */
private double[] guessAOmega(WeightedObservedPoint[] observations) {
    final double[] aOmega = new double[2];

    // initialize the sums for the linear model between the two integrals
    double sx2 = 0;
    double sy2 = 0;
    double sxy = 0;
    double sxz = 0;
    double syz = 0;

    double currentX = observations[0].getX();
    double currentY = observations[0].getY();
    double f2Integral = 0;
    double fPrime2Integral = 0;
    final double startX = currentX;
    for (int i = 1; i < observations.length; ++i) {
        // one step forward
        final double previousX = currentX;
        final double previousY = currentY;
        currentX = observations[i].getX();
        currentY = observations[i].getY();

        // update the integrals of f<sup>2</sup> and f'<sup>2</sup>
        // considering a linear model for f (and therefore constant f')
        final double dx = currentX - previousX;
        final double dy = currentY - previousY;
        final double f2StepIntegral =
            dx * (previousY * previousY + previousY * currentY + currentY * currentY) / 3;
        final double fPrime2StepIntegral = dy * dy / dx;

        final double x = currentX - startX;
        f2Integral += f2StepIntegral;
        fPrime2Integral += fPrime2StepIntegral;

        sx2 += x * x;
        sy2 += f2Integral * f2Integral;
        sxy += x * f2Integral;
        sxz += x * fPrime2Integral;
        syz += f2Integral * fPrime2Integral;
    }

    // compute the amplitude and pulsation coefficients
    double c1 = sy2 * sxz - sxy * syz;
    double c2 = sxy * sxz - sx2 * syz;
    double c3 = sx2 * sy2 - sxy * sxy;
    if ((c1 / c2 < 0) || (c2 / c3 < 0)) {
        final int last = observations.length - 1;
        // Range of the observations, assuming that the
        // observations are sorted.
        final double xRange = observations[last].getX() - observations[0].getX();
        if (xRange == 0) {
            throw new ZeroException();
        }
        aOmega[1] = 2 * Math.PI / xRange;

        double yMin = Double.POSITIVE_INFINITY;
        double yMax = Double.NEGATIVE_INFINITY;
        for (int i = 1; i < observations.length; ++i) {
            final double y = observations[i].getY();
            if (y < yMin) {
                yMin = y;
            }
            if (y > yMax) {
                yMax = y;
            }
        }
        aOmega[0] = 0.5 * (yMax - yMin);
    } else {
        if (c2 == 0) {
            // In some ill-conditioned cases (cf. MATH-844), the guesser
            // procedure cannot produce sensible results.
            throw new MathIllegalStateException(LocalizedFormats.ZERO_DENOMINATOR);
        }

        aOmega[0] = FastMath.sqrt(c1 / c2);
        aOmega[1] = FastMath.sqrt(c2 / c3);
    }

    return aOmega;
}
 

/**
 * <p>
 * Divide the value of this fraction by another, returning the result in
 * reduced form.
 * </p>
 *
 * @param fraction Fraction to divide by, must not be {@code null}.
 * @return a {@link BigFraction} instance with the resulting values.
 * @throws NullArgumentException if the {@code fraction} is {@code null}.
 * @throws MathArithmeticException if the fraction to divide by is zero
 */
public BigFraction divide(final BigFraction fraction) {
    if (fraction == null) {
        throw new NullArgumentException(LocalizedFormats.FRACTION);
    }
    if (BigInteger.ZERO.equals(fraction.numerator)) {
        throw new MathArithmeticException(LocalizedFormats.ZERO_DENOMINATOR);
    }

    return multiply(fraction.reciprocal());
}
 

/**
 * <p>
 * Divide the value of this fraction by the passed {@code BigInteger},
 * ie {@code this * 1 / bg}, returning the result in reduced form.
 * </p>
 *
 * @param bg the {@code BigInteger} to divide by, must not be {@code null}
 * @return a {@link BigFraction} instance with the resulting values
 * @throws NullArgumentException if the {@code BigInteger} is {@code null}
 * @throws MathArithmeticException if the fraction to divide by is zero
 */
public BigFraction divide(final BigInteger bg) {
    if (bg == null) {
        throw new NullArgumentException(LocalizedFormats.FRACTION);
    }
    if (BigInteger.ZERO.equals(bg)) {
        throw new MathArithmeticException(LocalizedFormats.ZERO_DENOMINATOR);
    }
    return new BigFraction(numerator, denominator.multiply(bg));
}
 

/**
 * <p>
 * Divide the value of this fraction by another, returning the result in
 * reduced form.
 * </p>
 *
 * @param fraction Fraction to divide by, must not be {@code null}.
 * @return a {@link BigFraction} instance with the resulting values.
 * @throws NullArgumentException if the {@code fraction} is {@code null}.
 * @throws MathArithmeticException if the fraction to divide by is zero
 */
public BigFraction divide(final BigFraction fraction) {
    if (fraction == null) {
        throw new NullArgumentException(LocalizedFormats.FRACTION);
    }
    if (BigInteger.ZERO.equals(fraction.numerator)) {
        throw new MathArithmeticException(LocalizedFormats.ZERO_DENOMINATOR);
    }

    return multiply(fraction.reciprocal());
}
 

/**
 * <p>
 * Divide the value of this fraction by the passed {@code BigInteger},
 * ie {@code this * 1 / bg}, returning the result in reduced form.
 * </p>
 *
 * @param bg the {@code BigInteger} to divide by, must not be {@code null}
 * @return a {@link BigFraction} instance with the resulting values
 * @throws NullArgumentException if the {@code BigInteger} is {@code null}
 * @throws MathArithmeticException if the fraction to divide by is zero
 */
public BigFraction divide(final BigInteger bg) {
    if (bg == null) {
        throw new NullArgumentException(LocalizedFormats.FRACTION);
    }
    if (BigInteger.ZERO.equals(bg)) {
        throw new MathArithmeticException(LocalizedFormats.ZERO_DENOMINATOR);
    }
    return new BigFraction(numerator, denominator.multiply(bg));
}
 

/**
 * <p>
 * Divide the value of this fraction by the passed {@code BigInteger},
 * ie {@code this * 1 / bg}, returning the result in reduced form.
 * </p>
 *
 * @param bg the {@code BigInteger} to divide by, must not be {@code null}
 * @return a {@link BigFraction} instance with the resulting values
 * @throws NullArgumentException if the {@code BigInteger} is {@code null}
 * @throws MathArithmeticException if the fraction to divide by is zero
 */
public BigFraction divide(final BigInteger bg) {
    if (bg == null) {
        throw new NullArgumentException(LocalizedFormats.FRACTION);
    }
    if (BigInteger.ZERO.equals(bg)) {
        throw new MathArithmeticException(LocalizedFormats.ZERO_DENOMINATOR);
    }
    return new BigFraction(numerator, denominator.multiply(bg));
}