You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by tn...@apache.org on 2015/02/16 23:40:07 UTC
[37/82] [partial] [math] Update for next development iteration:
commons-math4
http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/complex/ComplexFormat.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/complex/ComplexFormat.java b/src/main/java/org/apache/commons/math3/complex/ComplexFormat.java
deleted file mode 100644
index affb638..0000000
--- a/src/main/java/org/apache/commons/math3/complex/ComplexFormat.java
+++ /dev/null
@@ -1,427 +0,0 @@
-/*
- * Licensed to the Apache Software Foundation (ASF) under one or more
- * contributor license agreements. See the NOTICE file distributed with
- * this work for additional information regarding copyright ownership.
- * The ASF licenses this file to You under the Apache License, Version 2.0
- * (the "License"); you may not use this file except in compliance with
- * the License. You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-package org.apache.commons.math3.complex;
-
-import java.text.FieldPosition;
-import java.text.NumberFormat;
-import java.text.ParsePosition;
-import java.util.Locale;
-
-import org.apache.commons.math3.exception.MathIllegalArgumentException;
-import org.apache.commons.math3.exception.MathParseException;
-import org.apache.commons.math3.exception.NoDataException;
-import org.apache.commons.math3.exception.NullArgumentException;
-import org.apache.commons.math3.exception.util.LocalizedFormats;
-import org.apache.commons.math3.util.CompositeFormat;
-
-/**
- * Formats a Complex number in cartesian format "Re(c) + Im(c)i". 'i' can
- * be replaced with 'j' (or anything else), and the number format for both real
- * and imaginary parts can be configured.
- *
- */
-public class ComplexFormat {
-
- /** The default imaginary character. */
- private static final String DEFAULT_IMAGINARY_CHARACTER = "i";
- /** The notation used to signify the imaginary part of the complex number. */
- private final String imaginaryCharacter;
- /** The format used for the imaginary part. */
- private final NumberFormat imaginaryFormat;
- /** The format used for the real part. */
- private final NumberFormat realFormat;
-
- /**
- * Create an instance with the default imaginary character, 'i', and the
- * default number format for both real and imaginary parts.
- */
- public ComplexFormat() {
- this.imaginaryCharacter = DEFAULT_IMAGINARY_CHARACTER;
- this.imaginaryFormat = CompositeFormat.getDefaultNumberFormat();
- this.realFormat = imaginaryFormat;
- }
-
- /**
- * Create an instance with a custom number format for both real and
- * imaginary parts.
- * @param format the custom format for both real and imaginary parts.
- * @throws NullArgumentException if {@code realFormat} is {@code null}.
- */
- public ComplexFormat(NumberFormat format) throws NullArgumentException {
- if (format == null) {
- throw new NullArgumentException(LocalizedFormats.IMAGINARY_FORMAT);
- }
- this.imaginaryCharacter = DEFAULT_IMAGINARY_CHARACTER;
- this.imaginaryFormat = format;
- this.realFormat = format;
- }
-
- /**
- * Create an instance with a custom number format for the real part and a
- * custom number format for the imaginary part.
- * @param realFormat the custom format for the real part.
- * @param imaginaryFormat the custom format for the imaginary part.
- * @throws NullArgumentException if {@code imaginaryFormat} is {@code null}.
- * @throws NullArgumentException if {@code realFormat} is {@code null}.
- */
- public ComplexFormat(NumberFormat realFormat, NumberFormat imaginaryFormat)
- throws NullArgumentException {
- if (imaginaryFormat == null) {
- throw new NullArgumentException(LocalizedFormats.IMAGINARY_FORMAT);
- }
- if (realFormat == null) {
- throw new NullArgumentException(LocalizedFormats.REAL_FORMAT);
- }
-
- this.imaginaryCharacter = DEFAULT_IMAGINARY_CHARACTER;
- this.imaginaryFormat = imaginaryFormat;
- this.realFormat = realFormat;
- }
-
- /**
- * Create an instance with a custom imaginary character, and the default
- * number format for both real and imaginary parts.
- * @param imaginaryCharacter The custom imaginary character.
- * @throws NullArgumentException if {@code imaginaryCharacter} is
- * {@code null}.
- * @throws NoDataException if {@code imaginaryCharacter} is an
- * empty string.
- */
- public ComplexFormat(String imaginaryCharacter)
- throws NullArgumentException, NoDataException {
- this(imaginaryCharacter, CompositeFormat.getDefaultNumberFormat());
- }
-
- /**
- * Create an instance with a custom imaginary character, and a custom number
- * format for both real and imaginary parts.
- * @param imaginaryCharacter The custom imaginary character.
- * @param format the custom format for both real and imaginary parts.
- * @throws NullArgumentException if {@code imaginaryCharacter} is
- * {@code null}.
- * @throws NoDataException if {@code imaginaryCharacter} is an
- * empty string.
- * @throws NullArgumentException if {@code format} is {@code null}.
- */
- public ComplexFormat(String imaginaryCharacter, NumberFormat format)
- throws NullArgumentException, NoDataException {
- this(imaginaryCharacter, format, format);
- }
-
- /**
- * Create an instance with a custom imaginary character, a custom number
- * format for the real part, and a custom number format for the imaginary
- * part.
- *
- * @param imaginaryCharacter The custom imaginary character.
- * @param realFormat the custom format for the real part.
- * @param imaginaryFormat the custom format for the imaginary part.
- * @throws NullArgumentException if {@code imaginaryCharacter} is
- * {@code null}.
- * @throws NoDataException if {@code imaginaryCharacter} is an
- * empty string.
- * @throws NullArgumentException if {@code imaginaryFormat} is {@code null}.
- * @throws NullArgumentException if {@code realFormat} is {@code null}.
- */
- public ComplexFormat(String imaginaryCharacter,
- NumberFormat realFormat,
- NumberFormat imaginaryFormat)
- throws NullArgumentException, NoDataException {
- if (imaginaryCharacter == null) {
- throw new NullArgumentException();
- }
- if (imaginaryCharacter.length() == 0) {
- throw new NoDataException();
- }
- if (imaginaryFormat == null) {
- throw new NullArgumentException(LocalizedFormats.IMAGINARY_FORMAT);
- }
- if (realFormat == null) {
- throw new NullArgumentException(LocalizedFormats.REAL_FORMAT);
- }
-
- this.imaginaryCharacter = imaginaryCharacter;
- this.imaginaryFormat = imaginaryFormat;
- this.realFormat = realFormat;
- }
-
- /**
- * Get the set of locales for which complex formats are available.
- * <p>This is the same set as the {@link NumberFormat} set.</p>
- * @return available complex format locales.
- */
- public static Locale[] getAvailableLocales() {
- return NumberFormat.getAvailableLocales();
- }
-
- /**
- * This method calls {@link #format(Object,StringBuffer,FieldPosition)}.
- *
- * @param c Complex object to format.
- * @return A formatted number in the form "Re(c) + Im(c)i".
- */
- public String format(Complex c) {
- return format(c, new StringBuffer(), new FieldPosition(0)).toString();
- }
-
- /**
- * This method calls {@link #format(Object,StringBuffer,FieldPosition)}.
- *
- * @param c Double object to format.
- * @return A formatted number.
- */
- public String format(Double c) {
- return format(new Complex(c, 0), new StringBuffer(), new FieldPosition(0)).toString();
- }
-
- /**
- * Formats a {@link Complex} object to produce a string.
- *
- * @param complex the object to format.
- * @param toAppendTo where the text is to be appended
- * @param pos On input: an alignment field, if desired. On output: the
- * offsets of the alignment field
- * @return the value passed in as toAppendTo.
- */
- public StringBuffer format(Complex complex, StringBuffer toAppendTo,
- FieldPosition pos) {
- pos.setBeginIndex(0);
- pos.setEndIndex(0);
-
- // format real
- double re = complex.getReal();
- CompositeFormat.formatDouble(re, getRealFormat(), toAppendTo, pos);
-
- // format sign and imaginary
- double im = complex.getImaginary();
- StringBuffer imAppendTo;
- if (im < 0.0) {
- toAppendTo.append(" - ");
- imAppendTo = formatImaginary(-im, new StringBuffer(), pos);
- toAppendTo.append(imAppendTo);
- toAppendTo.append(getImaginaryCharacter());
- } else if (im > 0.0 || Double.isNaN(im)) {
- toAppendTo.append(" + ");
- imAppendTo = formatImaginary(im, new StringBuffer(), pos);
- toAppendTo.append(imAppendTo);
- toAppendTo.append(getImaginaryCharacter());
- }
-
- return toAppendTo;
- }
-
- /**
- * Format the absolute value of the imaginary part.
- *
- * @param absIm Absolute value of the imaginary part of a complex number.
- * @param toAppendTo where the text is to be appended.
- * @param pos On input: an alignment field, if desired. On output: the
- * offsets of the alignment field.
- * @return the value passed in as toAppendTo.
- */
- private StringBuffer formatImaginary(double absIm,
- StringBuffer toAppendTo,
- FieldPosition pos) {
- pos.setBeginIndex(0);
- pos.setEndIndex(0);
-
- CompositeFormat.formatDouble(absIm, getImaginaryFormat(), toAppendTo, pos);
- if (toAppendTo.toString().equals("1")) {
- // Remove the character "1" if it is the only one.
- toAppendTo.setLength(0);
- }
-
- return toAppendTo;
- }
-
- /**
- * Formats a object to produce a string. {@code obj} must be either a
- * {@link Complex} object or a {@link Number} object. Any other type of
- * object will result in an {@link IllegalArgumentException} being thrown.
- *
- * @param obj the object to format.
- * @param toAppendTo where the text is to be appended
- * @param pos On input: an alignment field, if desired. On output: the
- * offsets of the alignment field
- * @return the value passed in as toAppendTo.
- * @see java.text.Format#format(java.lang.Object, java.lang.StringBuffer, java.text.FieldPosition)
- * @throws MathIllegalArgumentException is {@code obj} is not a valid type.
- */
- public StringBuffer format(Object obj, StringBuffer toAppendTo,
- FieldPosition pos)
- throws MathIllegalArgumentException {
-
- StringBuffer ret = null;
-
- if (obj instanceof Complex) {
- ret = format( (Complex)obj, toAppendTo, pos);
- } else if (obj instanceof Number) {
- ret = format(new Complex(((Number)obj).doubleValue(), 0.0),
- toAppendTo, pos);
- } else {
- throw new MathIllegalArgumentException(LocalizedFormats.CANNOT_FORMAT_INSTANCE_AS_COMPLEX,
- obj.getClass().getName());
- }
-
- return ret;
- }
-
- /**
- * Access the imaginaryCharacter.
- * @return the imaginaryCharacter.
- */
- public String getImaginaryCharacter() {
- return imaginaryCharacter;
- }
-
- /**
- * Access the imaginaryFormat.
- * @return the imaginaryFormat.
- */
- public NumberFormat getImaginaryFormat() {
- return imaginaryFormat;
- }
-
- /**
- * Returns the default complex format for the current locale.
- * @return the default complex format.
- */
- public static ComplexFormat getInstance() {
- return getInstance(Locale.getDefault());
- }
-
- /**
- * Returns the default complex format for the given locale.
- * @param locale the specific locale used by the format.
- * @return the complex format specific to the given locale.
- */
- public static ComplexFormat getInstance(Locale locale) {
- NumberFormat f = CompositeFormat.getDefaultNumberFormat(locale);
- return new ComplexFormat(f);
- }
-
- /**
- * Returns the default complex format for the given locale.
- * @param locale the specific locale used by the format.
- * @param imaginaryCharacter Imaginary character.
- * @return the complex format specific to the given locale.
- * @throws NullArgumentException if {@code imaginaryCharacter} is
- * {@code null}.
- * @throws NoDataException if {@code imaginaryCharacter} is an
- * empty string.
- */
- public static ComplexFormat getInstance(String imaginaryCharacter, Locale locale)
- throws NullArgumentException, NoDataException {
- NumberFormat f = CompositeFormat.getDefaultNumberFormat(locale);
- return new ComplexFormat(imaginaryCharacter, f);
- }
-
- /**
- * Access the realFormat.
- * @return the realFormat.
- */
- public NumberFormat getRealFormat() {
- return realFormat;
- }
-
- /**
- * Parses a string to produce a {@link Complex} object.
- *
- * @param source the string to parse.
- * @return the parsed {@link Complex} object.
- * @throws MathParseException if the beginning of the specified string
- * cannot be parsed.
- */
- public Complex parse(String source) throws MathParseException {
- ParsePosition parsePosition = new ParsePosition(0);
- Complex result = parse(source, parsePosition);
- if (parsePosition.getIndex() == 0) {
- throw new MathParseException(source,
- parsePosition.getErrorIndex(),
- Complex.class);
- }
- return result;
- }
-
- /**
- * Parses a string to produce a {@link Complex} object.
- *
- * @param source the string to parse
- * @param pos input/ouput parsing parameter.
- * @return the parsed {@link Complex} object.
- */
- public Complex parse(String source, ParsePosition pos) {
- int initialIndex = pos.getIndex();
-
- // parse whitespace
- CompositeFormat.parseAndIgnoreWhitespace(source, pos);
-
- // parse real
- Number re = CompositeFormat.parseNumber(source, getRealFormat(), pos);
- if (re == null) {
- // invalid real number
- // set index back to initial, error index should already be set
- pos.setIndex(initialIndex);
- return null;
- }
-
- // parse sign
- int startIndex = pos.getIndex();
- char c = CompositeFormat.parseNextCharacter(source, pos);
- int sign = 0;
- switch (c) {
- case 0 :
- // no sign
- // return real only complex number
- return new Complex(re.doubleValue(), 0.0);
- case '-' :
- sign = -1;
- break;
- case '+' :
- sign = 1;
- break;
- default :
- // invalid sign
- // set index back to initial, error index should be the last
- // character examined.
- pos.setIndex(initialIndex);
- pos.setErrorIndex(startIndex);
- return null;
- }
-
- // parse whitespace
- CompositeFormat.parseAndIgnoreWhitespace(source, pos);
-
- // parse imaginary
- Number im = CompositeFormat.parseNumber(source, getRealFormat(), pos);
- if (im == null) {
- // invalid imaginary number
- // set index back to initial, error index should already be set
- pos.setIndex(initialIndex);
- return null;
- }
-
- // parse imaginary character
- if (!CompositeFormat.parseFixedstring(source, getImaginaryCharacter(), pos)) {
- return null;
- }
-
- return new Complex(re.doubleValue(), im.doubleValue() * sign);
-
- }
-}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/complex/ComplexUtils.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/complex/ComplexUtils.java b/src/main/java/org/apache/commons/math3/complex/ComplexUtils.java
deleted file mode 100644
index 55db946..0000000
--- a/src/main/java/org/apache/commons/math3/complex/ComplexUtils.java
+++ /dev/null
@@ -1,86 +0,0 @@
-/*
- * Licensed to the Apache Software Foundation (ASF) under one or more
- * contributor license agreements. See the NOTICE file distributed with
- * this work for additional information regarding copyright ownership.
- * The ASF licenses this file to You under the Apache License, Version 2.0
- * (the "License"); you may not use this file except in compliance with
- * the License. You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-package org.apache.commons.math3.complex;
-
-import org.apache.commons.math3.exception.MathIllegalArgumentException;
-import org.apache.commons.math3.exception.util.LocalizedFormats;
-import org.apache.commons.math3.util.FastMath;
-
-/**
- * Static implementations of common
- * {@link org.apache.commons.math3.complex.Complex} utilities functions.
- *
- */
-public class ComplexUtils {
-
- /**
- * Default constructor.
- */
- private ComplexUtils() {}
-
- /**
- * Creates a complex number from the given polar representation.
- * <p>
- * The value returned is <code>r·e<sup>i·theta</sup></code>,
- * computed as <code>r·cos(theta) + r·sin(theta)i</code></p>
- * <p>
- * If either <code>r</code> or <code>theta</code> is NaN, or
- * <code>theta</code> is infinite, {@link Complex#NaN} is returned.</p>
- * <p>
- * If <code>r</code> is infinite and <code>theta</code> is finite,
- * infinite or NaN values may be returned in parts of the result, following
- * the rules for double arithmetic.<pre>
- * Examples:
- * <code>
- * polar2Complex(INFINITY, π/4) = INFINITY + INFINITY i
- * polar2Complex(INFINITY, 0) = INFINITY + NaN i
- * polar2Complex(INFINITY, -π/4) = INFINITY - INFINITY i
- * polar2Complex(INFINITY, 5π/4) = -INFINITY - INFINITY i </code></pre></p>
- *
- * @param r the modulus of the complex number to create
- * @param theta the argument of the complex number to create
- * @return <code>r·e<sup>i·theta</sup></code>
- * @throws MathIllegalArgumentException if {@code r} is negative.
- * @since 1.1
- */
- public static Complex polar2Complex(double r, double theta) throws MathIllegalArgumentException {
- if (r < 0) {
- throw new MathIllegalArgumentException(
- LocalizedFormats.NEGATIVE_COMPLEX_MODULE, r);
- }
- return new Complex(r * FastMath.cos(theta), r * FastMath.sin(theta));
- }
-
- /**
- * Convert an array of primitive doubles to an array of {@code Complex} objects.
- *
- * @param real Array of numbers to be converted to their {@code Complex}
- * equivalent.
- * @return an array of {@code Complex} objects.
- *
- * @since 3.1
- */
- public static Complex[] convertToComplex(double[] real) {
- final Complex c[] = new Complex[real.length];
- for (int i = 0; i < real.length; i++) {
- c[i] = new Complex(real[i], 0);
- }
-
- return c;
- }
-}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/complex/Quaternion.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/complex/Quaternion.java b/src/main/java/org/apache/commons/math3/complex/Quaternion.java
deleted file mode 100644
index e845596..0000000
--- a/src/main/java/org/apache/commons/math3/complex/Quaternion.java
+++ /dev/null
@@ -1,464 +0,0 @@
-/*
- * Licensed to the Apache Software Foundation (ASF) under one or more
- * contributor license agreements. See the NOTICE file distributed with
- * this work for additional information regarding copyright ownership.
- * The ASF licenses this file to You under the Apache License, Version 2.0
- * (the "License"); you may not use this file except in compliance with
- * the License. You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-package org.apache.commons.math3.complex;
-
-import java.io.Serializable;
-import org.apache.commons.math3.util.FastMath;
-import org.apache.commons.math3.util.MathUtils;
-import org.apache.commons.math3.util.Precision;
-import org.apache.commons.math3.exception.DimensionMismatchException;
-import org.apache.commons.math3.exception.ZeroException;
-import org.apache.commons.math3.exception.util.LocalizedFormats;
-
-/**
- * This class implements <a href="http://mathworld.wolfram.com/Quaternion.html">
- * quaternions</a> (Hamilton's hypercomplex numbers).
- * <br/>
- * Instance of this class are guaranteed to be immutable.
- *
- * @since 3.1
- */
-public final class Quaternion implements Serializable {
- /** Identity quaternion. */
- public static final Quaternion IDENTITY = new Quaternion(1, 0, 0, 0);
- /** Zero quaternion. */
- public static final Quaternion ZERO = new Quaternion(0, 0, 0, 0);
- /** i */
- public static final Quaternion I = new Quaternion(0, 1, 0, 0);
- /** j */
- public static final Quaternion J = new Quaternion(0, 0, 1, 0);
- /** k */
- public static final Quaternion K = new Quaternion(0, 0, 0, 1);
-
- /** Serializable version identifier. */
- private static final long serialVersionUID = 20092012L;
-
- /** First component (scalar part). */
- private final double q0;
- /** Second component (first vector part). */
- private final double q1;
- /** Third component (second vector part). */
- private final double q2;
- /** Fourth component (third vector part). */
- private final double q3;
-
- /**
- * Builds a quaternion from its components.
- *
- * @param a Scalar component.
- * @param b First vector component.
- * @param c Second vector component.
- * @param d Third vector component.
- */
- public Quaternion(final double a,
- final double b,
- final double c,
- final double d) {
- this.q0 = a;
- this.q1 = b;
- this.q2 = c;
- this.q3 = d;
- }
-
- /**
- * Builds a quaternion from scalar and vector parts.
- *
- * @param scalar Scalar part of the quaternion.
- * @param v Components of the vector part of the quaternion.
- *
- * @throws DimensionMismatchException if the array length is not 3.
- */
- public Quaternion(final double scalar,
- final double[] v)
- throws DimensionMismatchException {
- if (v.length != 3) {
- throw new DimensionMismatchException(v.length, 3);
- }
- this.q0 = scalar;
- this.q1 = v[0];
- this.q2 = v[1];
- this.q3 = v[2];
- }
-
- /**
- * Builds a pure quaternion from a vector (assuming that the scalar
- * part is zero).
- *
- * @param v Components of the vector part of the pure quaternion.
- */
- public Quaternion(final double[] v) {
- this(0, v);
- }
-
- /**
- * Returns the conjugate quaternion of the instance.
- *
- * @return the conjugate quaternion
- */
- public Quaternion getConjugate() {
- return new Quaternion(q0, -q1, -q2, -q3);
- }
-
- /**
- * Returns the Hamilton product of two quaternions.
- *
- * @param q1 First quaternion.
- * @param q2 Second quaternion.
- * @return the product {@code q1} and {@code q2}, in that order.
- */
- public static Quaternion multiply(final Quaternion q1, final Quaternion q2) {
- // Components of the first quaternion.
- final double q1a = q1.getQ0();
- final double q1b = q1.getQ1();
- final double q1c = q1.getQ2();
- final double q1d = q1.getQ3();
-
- // Components of the second quaternion.
- final double q2a = q2.getQ0();
- final double q2b = q2.getQ1();
- final double q2c = q2.getQ2();
- final double q2d = q2.getQ3();
-
- // Components of the product.
- final double w = q1a * q2a - q1b * q2b - q1c * q2c - q1d * q2d;
- final double x = q1a * q2b + q1b * q2a + q1c * q2d - q1d * q2c;
- final double y = q1a * q2c - q1b * q2d + q1c * q2a + q1d * q2b;
- final double z = q1a * q2d + q1b * q2c - q1c * q2b + q1d * q2a;
-
- return new Quaternion(w, x, y, z);
- }
-
- /**
- * Returns the Hamilton product of the instance by a quaternion.
- *
- * @param q Quaternion.
- * @return the product of this instance with {@code q}, in that order.
- */
- public Quaternion multiply(final Quaternion q) {
- return multiply(this, q);
- }
-
- /**
- * Computes the sum of two quaternions.
- *
- * @param q1 Quaternion.
- * @param q2 Quaternion.
- * @return the sum of {@code q1} and {@code q2}.
- */
- public static Quaternion add(final Quaternion q1,
- final Quaternion q2) {
- return new Quaternion(q1.getQ0() + q2.getQ0(),
- q1.getQ1() + q2.getQ1(),
- q1.getQ2() + q2.getQ2(),
- q1.getQ3() + q2.getQ3());
- }
-
- /**
- * Computes the sum of the instance and another quaternion.
- *
- * @param q Quaternion.
- * @return the sum of this instance and {@code q}
- */
- public Quaternion add(final Quaternion q) {
- return add(this, q);
- }
-
- /**
- * Subtracts two quaternions.
- *
- * @param q1 First Quaternion.
- * @param q2 Second quaternion.
- * @return the difference between {@code q1} and {@code q2}.
- */
- public static Quaternion subtract(final Quaternion q1,
- final Quaternion q2) {
- return new Quaternion(q1.getQ0() - q2.getQ0(),
- q1.getQ1() - q2.getQ1(),
- q1.getQ2() - q2.getQ2(),
- q1.getQ3() - q2.getQ3());
- }
-
- /**
- * Subtracts a quaternion from the instance.
- *
- * @param q Quaternion.
- * @return the difference between this instance and {@code q}.
- */
- public Quaternion subtract(final Quaternion q) {
- return subtract(this, q);
- }
-
- /**
- * Computes the dot-product of two quaternions.
- *
- * @param q1 Quaternion.
- * @param q2 Quaternion.
- * @return the dot product of {@code q1} and {@code q2}.
- */
- public static double dotProduct(final Quaternion q1,
- final Quaternion q2) {
- return q1.getQ0() * q2.getQ0() +
- q1.getQ1() * q2.getQ1() +
- q1.getQ2() * q2.getQ2() +
- q1.getQ3() * q2.getQ3();
- }
-
- /**
- * Computes the dot-product of the instance by a quaternion.
- *
- * @param q Quaternion.
- * @return the dot product of this instance and {@code q}.
- */
- public double dotProduct(final Quaternion q) {
- return dotProduct(this, q);
- }
-
- /**
- * Computes the norm of the quaternion.
- *
- * @return the norm.
- */
- public double getNorm() {
- return FastMath.sqrt(q0 * q0 +
- q1 * q1 +
- q2 * q2 +
- q3 * q3);
- }
-
- /**
- * Computes the normalized quaternion (the versor of the instance).
- * The norm of the quaternion must not be zero.
- *
- * @return a normalized quaternion.
- * @throws ZeroException if the norm of the quaternion is zero.
- */
- public Quaternion normalize() {
- final double norm = getNorm();
-
- if (norm < Precision.SAFE_MIN) {
- throw new ZeroException(LocalizedFormats.NORM, norm);
- }
-
- return new Quaternion(q0 / norm,
- q1 / norm,
- q2 / norm,
- q3 / norm);
- }
-
- /**
- * {@inheritDoc}
- */
- @Override
- public boolean equals(Object other) {
- if (this == other) {
- return true;
- }
- if (other instanceof Quaternion) {
- final Quaternion q = (Quaternion) other;
- return q0 == q.getQ0() &&
- q1 == q.getQ1() &&
- q2 == q.getQ2() &&
- q3 == q.getQ3();
- }
-
- return false;
- }
-
- /**
- * {@inheritDoc}
- */
- @Override
- public int hashCode() {
- // "Effective Java" (second edition, p. 47).
- int result = 17;
- for (double comp : new double[] { q0, q1, q2, q3 }) {
- final int c = MathUtils.hash(comp);
- result = 31 * result + c;
- }
- return result;
- }
-
- /**
- * Checks whether this instance is equal to another quaternion
- * within a given tolerance.
- *
- * @param q Quaternion with which to compare the current quaternion.
- * @param eps Tolerance.
- * @return {@code true} if the each of the components are equal
- * within the allowed absolute error.
- */
- public boolean equals(final Quaternion q,
- final double eps) {
- return Precision.equals(q0, q.getQ0(), eps) &&
- Precision.equals(q1, q.getQ1(), eps) &&
- Precision.equals(q2, q.getQ2(), eps) &&
- Precision.equals(q3, q.getQ3(), eps);
- }
-
- /**
- * Checks whether the instance is a unit quaternion within a given
- * tolerance.
- *
- * @param eps Tolerance (absolute error).
- * @return {@code true} if the norm is 1 within the given tolerance,
- * {@code false} otherwise
- */
- public boolean isUnitQuaternion(double eps) {
- return Precision.equals(getNorm(), 1d, eps);
- }
-
- /**
- * Checks whether the instance is a pure quaternion within a given
- * tolerance.
- *
- * @param eps Tolerance (absolute error).
- * @return {@code true} if the scalar part of the quaternion is zero.
- */
- public boolean isPureQuaternion(double eps) {
- return FastMath.abs(getQ0()) <= eps;
- }
-
- /**
- * Returns the polar form of the quaternion.
- *
- * @return the unit quaternion with positive scalar part.
- */
- public Quaternion getPositivePolarForm() {
- if (getQ0() < 0) {
- final Quaternion unitQ = normalize();
- // The quaternion of rotation (normalized quaternion) q and -q
- // are equivalent (i.e. represent the same rotation).
- return new Quaternion(-unitQ.getQ0(),
- -unitQ.getQ1(),
- -unitQ.getQ2(),
- -unitQ.getQ3());
- } else {
- return this.normalize();
- }
- }
-
- /**
- * Returns the inverse of this instance.
- * The norm of the quaternion must not be zero.
- *
- * @return the inverse.
- * @throws ZeroException if the norm (squared) of the quaternion is zero.
- */
- public Quaternion getInverse() {
- final double squareNorm = q0 * q0 + q1 * q1 + q2 * q2 + q3 * q3;
- if (squareNorm < Precision.SAFE_MIN) {
- throw new ZeroException(LocalizedFormats.NORM, squareNorm);
- }
-
- return new Quaternion(q0 / squareNorm,
- -q1 / squareNorm,
- -q2 / squareNorm,
- -q3 / squareNorm);
- }
-
- /**
- * Gets the first component of the quaternion (scalar part).
- *
- * @return the scalar part.
- */
- public double getQ0() {
- return q0;
- }
-
- /**
- * Gets the second component of the quaternion (first component
- * of the vector part).
- *
- * @return the first component of the vector part.
- */
- public double getQ1() {
- return q1;
- }
-
- /**
- * Gets the third component of the quaternion (second component
- * of the vector part).
- *
- * @return the second component of the vector part.
- */
- public double getQ2() {
- return q2;
- }
-
- /**
- * Gets the fourth component of the quaternion (third component
- * of the vector part).
- *
- * @return the third component of the vector part.
- */
- public double getQ3() {
- return q3;
- }
-
- /**
- * Gets the scalar part of the quaternion.
- *
- * @return the scalar part.
- * @see #getQ0()
- */
- public double getScalarPart() {
- return getQ0();
- }
-
- /**
- * Gets the three components of the vector part of the quaternion.
- *
- * @return the vector part.
- * @see #getQ1()
- * @see #getQ2()
- * @see #getQ3()
- */
- public double[] getVectorPart() {
- return new double[] { getQ1(), getQ2(), getQ3() };
- }
-
- /**
- * Multiplies the instance by a scalar.
- *
- * @param alpha Scalar factor.
- * @return a scaled quaternion.
- */
- public Quaternion multiply(final double alpha) {
- return new Quaternion(alpha * q0,
- alpha * q1,
- alpha * q2,
- alpha * q3);
- }
-
- /**
- * {@inheritDoc}
- */
- @Override
- public String toString() {
- final String sp = " ";
- final StringBuilder s = new StringBuilder();
- s.append("[")
- .append(q0).append(sp)
- .append(q1).append(sp)
- .append(q2).append(sp)
- .append(q3)
- .append("]");
-
- return s.toString();
- }
-}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/complex/RootsOfUnity.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/complex/RootsOfUnity.java b/src/main/java/org/apache/commons/math3/complex/RootsOfUnity.java
deleted file mode 100644
index 4e63835..0000000
--- a/src/main/java/org/apache/commons/math3/complex/RootsOfUnity.java
+++ /dev/null
@@ -1,218 +0,0 @@
-/*
- * Licensed to the Apache Software Foundation (ASF) under one or more
- * contributor license agreements. See the NOTICE file distributed with
- * this work for additional information regarding copyright ownership.
- * The ASF licenses this file to You under the Apache License, Version 2.0
- * (the "License"); you may not use this file except in compliance with
- * the License. You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-package org.apache.commons.math3.complex;
-
-import java.io.Serializable;
-
-import org.apache.commons.math3.exception.MathIllegalArgumentException;
-import org.apache.commons.math3.exception.MathIllegalStateException;
-import org.apache.commons.math3.exception.OutOfRangeException;
-import org.apache.commons.math3.exception.ZeroException;
-import org.apache.commons.math3.exception.util.LocalizedFormats;
-import org.apache.commons.math3.util.FastMath;
-
-/**
- * A helper class for the computation and caching of the {@code n}-th roots of
- * unity.
- *
- * @since 3.0
- */
-public class RootsOfUnity implements Serializable {
-
- /** Serializable version id. */
- private static final long serialVersionUID = 20120201L;
-
- /** Number of roots of unity. */
- private int omegaCount;
-
- /** Real part of the roots. */
- private double[] omegaReal;
-
- /**
- * Imaginary part of the {@code n}-th roots of unity, for positive values
- * of {@code n}. In this array, the roots are stored in counter-clockwise
- * order.
- */
- private double[] omegaImaginaryCounterClockwise;
-
- /**
- * Imaginary part of the {@code n}-th roots of unity, for negative values
- * of {@code n}. In this array, the roots are stored in clockwise order.
- */
- private double[] omegaImaginaryClockwise;
-
- /**
- * {@code true} if {@link #computeRoots(int)} was called with a positive
- * value of its argument {@code n}. In this case, counter-clockwise ordering
- * of the roots of unity should be used.
- */
- private boolean isCounterClockWise;
-
- /**
- * Build an engine for computing the {@code n}-th roots of unity.
- */
- public RootsOfUnity() {
-
- omegaCount = 0;
- omegaReal = null;
- omegaImaginaryCounterClockwise = null;
- omegaImaginaryClockwise = null;
- isCounterClockWise = true;
- }
-
- /**
- * Returns {@code true} if {@link #computeRoots(int)} was called with a
- * positive value of its argument {@code n}. If {@code true}, then
- * counter-clockwise ordering of the roots of unity should be used.
- *
- * @return {@code true} if the roots of unity are stored in
- * counter-clockwise order
- * @throws MathIllegalStateException if no roots of unity have been computed
- * yet
- */
- public synchronized boolean isCounterClockWise()
- throws MathIllegalStateException {
-
- if (omegaCount == 0) {
- throw new MathIllegalStateException(
- LocalizedFormats.ROOTS_OF_UNITY_NOT_COMPUTED_YET);
- }
- return isCounterClockWise;
- }
-
- /**
- * <p>
- * Computes the {@code n}-th roots of unity. The roots are stored in
- * {@code omega[]}, such that {@code omega[k] = w ^ k}, where
- * {@code k = 0, ..., n - 1}, {@code w = exp(2 * pi * i / n)} and
- * {@code i = sqrt(-1)}.
- * </p>
- * <p>
- * Note that {@code n} can be positive of negative
- * </p>
- * <ul>
- * <li>{@code abs(n)} is always the number of roots of unity.</li>
- * <li>If {@code n > 0}, then the roots are stored in counter-clockwise order.</li>
- * <li>If {@code n < 0}, then the roots are stored in clockwise order.</p>
- * </ul>
- *
- * @param n the (signed) number of roots of unity to be computed
- * @throws ZeroException if {@code n = 0}
- */
- public synchronized void computeRoots(int n) throws ZeroException {
-
- if (n == 0) {
- throw new ZeroException(
- LocalizedFormats.CANNOT_COMPUTE_0TH_ROOT_OF_UNITY);
- }
-
- isCounterClockWise = n > 0;
-
- // avoid repetitive calculations
- final int absN = FastMath.abs(n);
-
- if (absN == omegaCount) {
- return;
- }
-
- // calculate everything from scratch
- final double t = 2.0 * FastMath.PI / absN;
- final double cosT = FastMath.cos(t);
- final double sinT = FastMath.sin(t);
- omegaReal = new double[absN];
- omegaImaginaryCounterClockwise = new double[absN];
- omegaImaginaryClockwise = new double[absN];
- omegaReal[0] = 1.0;
- omegaImaginaryCounterClockwise[0] = 0.0;
- omegaImaginaryClockwise[0] = 0.0;
- for (int i = 1; i < absN; i++) {
- omegaReal[i] = omegaReal[i - 1] * cosT -
- omegaImaginaryCounterClockwise[i - 1] * sinT;
- omegaImaginaryCounterClockwise[i] = omegaReal[i - 1] * sinT +
- omegaImaginaryCounterClockwise[i - 1] * cosT;
- omegaImaginaryClockwise[i] = -omegaImaginaryCounterClockwise[i];
- }
- omegaCount = absN;
- }
-
- /**
- * Get the real part of the {@code k}-th {@code n}-th root of unity.
- *
- * @param k index of the {@code n}-th root of unity
- * @return real part of the {@code k}-th {@code n}-th root of unity
- * @throws MathIllegalStateException if no roots of unity have been
- * computed yet
- * @throws MathIllegalArgumentException if {@code k} is out of range
- */
- public synchronized double getReal(int k)
- throws MathIllegalStateException, MathIllegalArgumentException {
-
- if (omegaCount == 0) {
- throw new MathIllegalStateException(
- LocalizedFormats.ROOTS_OF_UNITY_NOT_COMPUTED_YET);
- }
- if ((k < 0) || (k >= omegaCount)) {
- throw new OutOfRangeException(
- LocalizedFormats.OUT_OF_RANGE_ROOT_OF_UNITY_INDEX,
- Integer.valueOf(k),
- Integer.valueOf(0),
- Integer.valueOf(omegaCount - 1));
- }
-
- return omegaReal[k];
- }
-
- /**
- * Get the imaginary part of the {@code k}-th {@code n}-th root of unity.
- *
- * @param k index of the {@code n}-th root of unity
- * @return imaginary part of the {@code k}-th {@code n}-th root of unity
- * @throws MathIllegalStateException if no roots of unity have been
- * computed yet
- * @throws OutOfRangeException if {@code k} is out of range
- */
- public synchronized double getImaginary(int k)
- throws MathIllegalStateException, OutOfRangeException {
-
- if (omegaCount == 0) {
- throw new MathIllegalStateException(
- LocalizedFormats.ROOTS_OF_UNITY_NOT_COMPUTED_YET);
- }
- if ((k < 0) || (k >= omegaCount)) {
- throw new OutOfRangeException(
- LocalizedFormats.OUT_OF_RANGE_ROOT_OF_UNITY_INDEX,
- Integer.valueOf(k),
- Integer.valueOf(0),
- Integer.valueOf(omegaCount - 1));
- }
-
- return isCounterClockWise ? omegaImaginaryCounterClockwise[k] :
- omegaImaginaryClockwise[k];
- }
-
- /**
- * Returns the number of roots of unity currently stored. If
- * {@link #computeRoots(int)} was called with {@code n}, then this method
- * returns {@code abs(n)}. If no roots of unity have been computed yet, this
- * method returns 0.
- *
- * @return the number of roots of unity currently stored
- */
- public synchronized int getNumberOfRoots() {
- return omegaCount;
- }
-}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/complex/package-info.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/complex/package-info.java b/src/main/java/org/apache/commons/math3/complex/package-info.java
deleted file mode 100644
index 818806d..0000000
--- a/src/main/java/org/apache/commons/math3/complex/package-info.java
+++ /dev/null
@@ -1,23 +0,0 @@
-/*
- * 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.
- */
-/**
- *
- * Complex number type and implementations of complex transcendental
- * functions.
- *
- */
-package org.apache.commons.math3.complex;
http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/dfp/BracketingNthOrderBrentSolverDFP.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/dfp/BracketingNthOrderBrentSolverDFP.java b/src/main/java/org/apache/commons/math3/dfp/BracketingNthOrderBrentSolverDFP.java
deleted file mode 100644
index c0328b6..0000000
--- a/src/main/java/org/apache/commons/math3/dfp/BracketingNthOrderBrentSolverDFP.java
+++ /dev/null
@@ -1,436 +0,0 @@
-/*
- * Licensed to the Apache Software Foundation (ASF) under one or more
- * contributor license agreements. See the NOTICE file distributed with
- * this work for additional information regarding copyright ownership.
- * The ASF licenses this file to You under the Apache License, Version 2.0
- * (the "License"); you may not use this file except in compliance with
- * the License. You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-package org.apache.commons.math3.dfp;
-
-
-import org.apache.commons.math3.analysis.solvers.AllowedSolution;
-import org.apache.commons.math3.exception.MathInternalError;
-import org.apache.commons.math3.exception.NoBracketingException;
-import org.apache.commons.math3.exception.NullArgumentException;
-import org.apache.commons.math3.exception.NumberIsTooSmallException;
-import org.apache.commons.math3.util.Incrementor;
-import org.apache.commons.math3.util.MathUtils;
-
-/**
- * This class implements a modification of the <a
- * href="http://mathworld.wolfram.com/BrentsMethod.html"> Brent algorithm</a>.
- * <p>
- * The changes with respect to the original Brent algorithm are:
- * <ul>
- * <li>the returned value is chosen in the current interval according
- * to user specified {@link AllowedSolution},</li>
- * <li>the maximal order for the invert polynomial root search is
- * user-specified instead of being invert quadratic only</li>
- * </ul>
- * </p>
- * The given interval must bracket the root.
- *
- */
-public class BracketingNthOrderBrentSolverDFP {
-
- /** Maximal aging triggering an attempt to balance the bracketing interval. */
- private static final int MAXIMAL_AGING = 2;
-
- /** Maximal order. */
- private final int maximalOrder;
-
- /** Function value accuracy. */
- private final Dfp functionValueAccuracy;
-
- /** Absolute accuracy. */
- private final Dfp absoluteAccuracy;
-
- /** Relative accuracy. */
- private final Dfp relativeAccuracy;
-
- /** Evaluations counter. */
- private final Incrementor evaluations = new Incrementor();
-
- /**
- * Construct a solver.
- *
- * @param relativeAccuracy Relative accuracy.
- * @param absoluteAccuracy Absolute accuracy.
- * @param functionValueAccuracy Function value accuracy.
- * @param maximalOrder maximal order.
- * @exception NumberIsTooSmallException if maximal order is lower than 2
- */
- public BracketingNthOrderBrentSolverDFP(final Dfp relativeAccuracy,
- final Dfp absoluteAccuracy,
- final Dfp functionValueAccuracy,
- final int maximalOrder)
- throws NumberIsTooSmallException {
- if (maximalOrder < 2) {
- throw new NumberIsTooSmallException(maximalOrder, 2, true);
- }
- this.maximalOrder = maximalOrder;
- this.absoluteAccuracy = absoluteAccuracy;
- this.relativeAccuracy = relativeAccuracy;
- this.functionValueAccuracy = functionValueAccuracy;
- }
-
- /** Get the maximal order.
- * @return maximal order
- */
- public int getMaximalOrder() {
- return maximalOrder;
- }
-
- /**
- * Get the maximal number of function evaluations.
- *
- * @return the maximal number of function evaluations.
- */
- public int getMaxEvaluations() {
- return evaluations.getMaximalCount();
- }
-
- /**
- * Get the number of evaluations of the objective function.
- * The number of evaluations corresponds to the last call to the
- * {@code optimize} method. It is 0 if the method has not been
- * called yet.
- *
- * @return the number of evaluations of the objective function.
- */
- public int getEvaluations() {
- return evaluations.getCount();
- }
-
- /**
- * Get the absolute accuracy.
- * @return absolute accuracy
- */
- public Dfp getAbsoluteAccuracy() {
- return absoluteAccuracy;
- }
-
- /**
- * Get the relative accuracy.
- * @return relative accuracy
- */
- public Dfp getRelativeAccuracy() {
- return relativeAccuracy;
- }
-
- /**
- * Get the function accuracy.
- * @return function accuracy
- */
- public Dfp getFunctionValueAccuracy() {
- return functionValueAccuracy;
- }
-
- /**
- * Solve for a zero in the given interval.
- * A solver may require that the interval brackets a single zero root.
- * Solvers that do require bracketing should be able to handle the case
- * where one of the endpoints is itself a root.
- *
- * @param maxEval Maximum number of evaluations.
- * @param f Function to solve.
- * @param min Lower bound for the interval.
- * @param max Upper bound for the interval.
- * @param allowedSolution The kind of solutions that the root-finding algorithm may
- * accept as solutions.
- * @return a value where the function is zero.
- * @exception NullArgumentException if f is null.
- * @exception NoBracketingException if root cannot be bracketed
- */
- public Dfp solve(final int maxEval, final UnivariateDfpFunction f,
- final Dfp min, final Dfp max, final AllowedSolution allowedSolution)
- throws NullArgumentException, NoBracketingException {
- return solve(maxEval, f, min, max, min.add(max).divide(2), allowedSolution);
- }
-
- /**
- * Solve for a zero in the given interval, start at {@code startValue}.
- * A solver may require that the interval brackets a single zero root.
- * Solvers that do require bracketing should be able to handle the case
- * where one of the endpoints is itself a root.
- *
- * @param maxEval Maximum number of evaluations.
- * @param f Function to solve.
- * @param min Lower bound for the interval.
- * @param max Upper bound for the interval.
- * @param startValue Start value to use.
- * @param allowedSolution The kind of solutions that the root-finding algorithm may
- * accept as solutions.
- * @return a value where the function is zero.
- * @exception NullArgumentException if f is null.
- * @exception NoBracketingException if root cannot be bracketed
- */
- public Dfp solve(final int maxEval, final UnivariateDfpFunction f,
- final Dfp min, final Dfp max, final Dfp startValue,
- final AllowedSolution allowedSolution)
- throws NullArgumentException, NoBracketingException {
-
- // Checks.
- MathUtils.checkNotNull(f);
-
- // Reset.
- evaluations.setMaximalCount(maxEval);
- evaluations.resetCount();
- Dfp zero = startValue.getZero();
- Dfp nan = zero.newInstance((byte) 1, Dfp.QNAN);
-
- // prepare arrays with the first points
- final Dfp[] x = new Dfp[maximalOrder + 1];
- final Dfp[] y = new Dfp[maximalOrder + 1];
- x[0] = min;
- x[1] = startValue;
- x[2] = max;
-
- // evaluate initial guess
- evaluations.incrementCount();
- y[1] = f.value(x[1]);
- if (y[1].isZero()) {
- // return the initial guess if it is a perfect root.
- return x[1];
- }
-
- // evaluate first endpoint
- evaluations.incrementCount();
- y[0] = f.value(x[0]);
- if (y[0].isZero()) {
- // return the first endpoint if it is a perfect root.
- return x[0];
- }
-
- int nbPoints;
- int signChangeIndex;
- if (y[0].multiply(y[1]).negativeOrNull()) {
-
- // reduce interval if it brackets the root
- nbPoints = 2;
- signChangeIndex = 1;
-
- } else {
-
- // evaluate second endpoint
- evaluations.incrementCount();
- y[2] = f.value(x[2]);
- if (y[2].isZero()) {
- // return the second endpoint if it is a perfect root.
- return x[2];
- }
-
- if (y[1].multiply(y[2]).negativeOrNull()) {
- // use all computed point as a start sampling array for solving
- nbPoints = 3;
- signChangeIndex = 2;
- } else {
- throw new NoBracketingException(x[0].toDouble(), x[2].toDouble(),
- y[0].toDouble(), y[2].toDouble());
- }
-
- }
-
- // prepare a work array for inverse polynomial interpolation
- final Dfp[] tmpX = new Dfp[x.length];
-
- // current tightest bracketing of the root
- Dfp xA = x[signChangeIndex - 1];
- Dfp yA = y[signChangeIndex - 1];
- Dfp absXA = xA.abs();
- Dfp absYA = yA.abs();
- int agingA = 0;
- Dfp xB = x[signChangeIndex];
- Dfp yB = y[signChangeIndex];
- Dfp absXB = xB.abs();
- Dfp absYB = yB.abs();
- int agingB = 0;
-
- // search loop
- while (true) {
-
- // check convergence of bracketing interval
- Dfp maxX = absXA.lessThan(absXB) ? absXB : absXA;
- Dfp maxY = absYA.lessThan(absYB) ? absYB : absYA;
- final Dfp xTol = absoluteAccuracy.add(relativeAccuracy.multiply(maxX));
- if (xB.subtract(xA).subtract(xTol).negativeOrNull() ||
- maxY.lessThan(functionValueAccuracy)) {
- switch (allowedSolution) {
- case ANY_SIDE :
- return absYA.lessThan(absYB) ? xA : xB;
- case LEFT_SIDE :
- return xA;
- case RIGHT_SIDE :
- return xB;
- case BELOW_SIDE :
- return yA.lessThan(zero) ? xA : xB;
- case ABOVE_SIDE :
- return yA.lessThan(zero) ? xB : xA;
- default :
- // this should never happen
- throw new MathInternalError(null);
- }
- }
-
- // target for the next evaluation point
- Dfp targetY;
- if (agingA >= MAXIMAL_AGING) {
- // we keep updating the high bracket, try to compensate this
- targetY = yB.divide(16).negate();
- } else if (agingB >= MAXIMAL_AGING) {
- // we keep updating the low bracket, try to compensate this
- targetY = yA.divide(16).negate();
- } else {
- // bracketing is balanced, try to find the root itself
- targetY = zero;
- }
-
- // make a few attempts to guess a root,
- Dfp nextX;
- int start = 0;
- int end = nbPoints;
- do {
-
- // guess a value for current target, using inverse polynomial interpolation
- System.arraycopy(x, start, tmpX, start, end - start);
- nextX = guessX(targetY, tmpX, y, start, end);
-
- if (!(nextX.greaterThan(xA) && nextX.lessThan(xB))) {
- // the guessed root is not strictly inside of the tightest bracketing interval
-
- // the guessed root is either not strictly inside the interval or it
- // is a NaN (which occurs when some sampling points share the same y)
- // we try again with a lower interpolation order
- if (signChangeIndex - start >= end - signChangeIndex) {
- // we have more points before the sign change, drop the lowest point
- ++start;
- } else {
- // we have more points after sign change, drop the highest point
- --end;
- }
-
- // we need to do one more attempt
- nextX = nan;
-
- }
-
- } while (nextX.isNaN() && (end - start > 1));
-
- if (nextX.isNaN()) {
- // fall back to bisection
- nextX = xA.add(xB.subtract(xA).divide(2));
- start = signChangeIndex - 1;
- end = signChangeIndex;
- }
-
- // evaluate the function at the guessed root
- evaluations.incrementCount();
- final Dfp nextY = f.value(nextX);
- if (nextY.isZero()) {
- // we have found an exact root, since it is not an approximation
- // we don't need to bother about the allowed solutions setting
- return nextX;
- }
-
- if ((nbPoints > 2) && (end - start != nbPoints)) {
-
- // we have been forced to ignore some points to keep bracketing,
- // they are probably too far from the root, drop them from now on
- nbPoints = end - start;
- System.arraycopy(x, start, x, 0, nbPoints);
- System.arraycopy(y, start, y, 0, nbPoints);
- signChangeIndex -= start;
-
- } else if (nbPoints == x.length) {
-
- // we have to drop one point in order to insert the new one
- nbPoints--;
-
- // keep the tightest bracketing interval as centered as possible
- if (signChangeIndex >= (x.length + 1) / 2) {
- // we drop the lowest point, we have to shift the arrays and the index
- System.arraycopy(x, 1, x, 0, nbPoints);
- System.arraycopy(y, 1, y, 0, nbPoints);
- --signChangeIndex;
- }
-
- }
-
- // insert the last computed point
- //(by construction, we know it lies inside the tightest bracketing interval)
- System.arraycopy(x, signChangeIndex, x, signChangeIndex + 1, nbPoints - signChangeIndex);
- x[signChangeIndex] = nextX;
- System.arraycopy(y, signChangeIndex, y, signChangeIndex + 1, nbPoints - signChangeIndex);
- y[signChangeIndex] = nextY;
- ++nbPoints;
-
- // update the bracketing interval
- if (nextY.multiply(yA).negativeOrNull()) {
- // the sign change occurs before the inserted point
- xB = nextX;
- yB = nextY;
- absYB = yB.abs();
- ++agingA;
- agingB = 0;
- } else {
- // the sign change occurs after the inserted point
- xA = nextX;
- yA = nextY;
- absYA = yA.abs();
- agingA = 0;
- ++agingB;
-
- // update the sign change index
- signChangeIndex++;
-
- }
-
- }
-
- }
-
- /** Guess an x value by n<sup>th</sup> order inverse polynomial interpolation.
- * <p>
- * The x value is guessed by evaluating polynomial Q(y) at y = targetY, where Q
- * is built such that for all considered points (x<sub>i</sub>, y<sub>i</sub>),
- * Q(y<sub>i</sub>) = x<sub>i</sub>.
- * </p>
- * @param targetY target value for y
- * @param x reference points abscissas for interpolation,
- * note that this array <em>is</em> modified during computation
- * @param y reference points ordinates for interpolation
- * @param start start index of the points to consider (inclusive)
- * @param end end index of the points to consider (exclusive)
- * @return guessed root (will be a NaN if two points share the same y)
- */
- private Dfp guessX(final Dfp targetY, final Dfp[] x, final Dfp[] y,
- final int start, final int end) {
-
- // compute Q Newton coefficients by divided differences
- for (int i = start; i < end - 1; ++i) {
- final int delta = i + 1 - start;
- for (int j = end - 1; j > i; --j) {
- x[j] = x[j].subtract(x[j-1]).divide(y[j].subtract(y[j - delta]));
- }
- }
-
- // evaluate Q(targetY)
- Dfp x0 = targetY.getZero();
- for (int j = end - 1; j >= start; --j) {
- x0 = x[j].add(x0.multiply(targetY.subtract(y[j])));
- }
-
- return x0;
-
- }
-
-}