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; - - } - -}
