| /* |
| * 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.sis.math; |
| |
| import java.util.Arrays; |
| import org.apache.sis.util.Static; |
| import org.apache.sis.util.ArraysExt; |
| import org.apache.sis.util.ArgumentChecks; |
| import org.apache.sis.util.resources.Errors; |
| import org.apache.sis.internal.util.DoubleDouble; |
| |
| import static java.lang.Math.abs; |
| import static java.lang.Math.sqrt; |
| import static java.lang.Float.intBitsToFloat; |
| import static java.lang.Float.floatToIntBits; |
| import static java.lang.Float.floatToRawIntBits; |
| import static java.lang.Double.longBitsToDouble; |
| import static java.lang.Double.doubleToLongBits; |
| import static java.lang.Double.doubleToRawLongBits; |
| import static org.apache.sis.internal.util.Numerics.SIGN_BIT_MASK; |
| import static org.apache.sis.internal.util.Numerics.SIGNIFICAND_SIZE; |
| |
| |
| /** |
| * Simple mathematical functions in addition to the ones provided in {@link Math}. |
| * Some methods in this class are very similar to the standard {@link Math} methods |
| * or could be implemented with straightforward formulas. |
| * However the methods in this class put an emphasis on: |
| * |
| * <ul> |
| * <li>Rounding errors: |
| * {@link #magnitude(double[]) magnitude}, |
| * {@link #pow10(double) pow10}.</li> |
| * <li>Distinguishing positive zero from negative zero: |
| * {@link #isPositive(double) isPositive}, |
| * {@link #isNegative(double) isNegative}, |
| * {@link #isSameSign(double, double) isSameSign}, |
| * {@link #xorSign(double, double) xorSign}.</li> |
| * <li>Distinguishing the different kinds of NaN numbers: |
| * {@link #toNanFloat(int) toNanFloat}, |
| * {@link #toNanOrdinal(float) toNanOrdinal}.</li> |
| * </ul> |
| * |
| * Some additional functions not found in {@code Math} are: |
| * {@link #atanh(double) atanh}, |
| * {@link #nextPrimeNumber(int) nextPrimeNumber}. |
| * |
| * @author Martin Desruisseaux (MPO, IRD, Geomatys) |
| * @author Johann Sorel (Geomatys) |
| * @version 1.0 |
| * |
| * @see DecimalFunctions |
| * @see org.apache.sis.util.Numbers |
| * |
| * @since 0.3 |
| * @module |
| */ |
| public final class MathFunctions extends Static { |
| /** |
| * The square root of 2, which is approximated by {@value}. |
| * |
| * @see Math#sqrt(double) |
| */ |
| public static final double SQRT_2 = 1.4142135623730951; |
| |
| /** |
| * The logarithm of 2 in base 10, which is approximated by {@value}. |
| * This constant is useful for converting a power of 2 to a power of 10 as below: |
| * |
| * {@preformat java |
| * double exp10 = exp2 * LOG10_2; |
| * } |
| * |
| * @see Math#log10(double) |
| * @see #getExponent(double) |
| * |
| * @since 0.4 |
| */ |
| public static final double LOG10_2 = 0.3010299956639812; |
| |
| /** |
| * The minimal ordinal value for {@code NaN} numbers created by {@link #toNanFloat(int)}. |
| * The current value is {@value}. |
| * |
| * @since 1.0 |
| */ |
| public static final int MIN_NAN_ORDINAL = -0x200000; |
| |
| /** |
| * The maximal ordinal value for {@code NaN} numbers created by {@link #toNanFloat(int)}. |
| * The current value is {@value}. |
| * |
| * @since 1.0 |
| */ |
| public static final int MAX_NAN_ORDINAL = 0x1FFFFF; |
| |
| /** |
| * The beginning of ranges of quiet NaN values. |
| * The range is selected in way to restrict ourselves to <cite>quiet</cite> NaN values. |
| * The following is an adaptation of evaluator's comments for bug #4471414 |
| * (http://developer.java.sun.com/developer/bugParade/bugs/4471414.html): |
| * |
| * <blockquote> |
| * There are actually two types of NaNs, signaling NaNs and quiet NaNs. Java doesn't support the features necessary |
| * to reliably distinguish the two. However, the relevant point is that copying a signaling NaN may (or may not, at |
| * the implementers discretion) yield a quiet NaN — a NaN with a different bit pattern (IEEE 754 6.2). Therefore, on |
| * IEEE 754 compliant platforms it may be impossible to find a signaling NaN stored in an array since a signaling NaN |
| * passed as an argument to binarySearch may get replaced by a quiet NaN. |
| * </blockquote> |
| * |
| * The relevant thresholds are: |
| * <ul> |
| * <li>{@code 7F800000}: positive infinity.</li> |
| * <li>{@code 7F800001}: first signaling NaN in the range of positive values.</li> |
| * <li>{@code 7FBFFFFF}: last signaling NaN.</li> |
| * <li>{@code 7FC00000}: first quiet NaN. Also the standard {@link Double#NaN} value.</li> |
| * <li>{@code 7FFFFFFF}: last quiet NaN.</li> |
| * <li>{@code FF800000}: negative infinity.</li> |
| * <li>{@code FF800001}: first signaling NaN in the range of negative values.</li> |
| * <li>{@code FFBFFFFF}: last signaling NaN.</li> |
| * <li>{@code FFC00000}: first quiet NaN in the range of negative values.</li> |
| * <li>{@code FFFFFFFF}: last quiet NaN.</li> |
| * </ul> |
| * |
| * @see #toNanFloat(int) |
| * @see #toNanOrdinal(float) |
| */ |
| static final int POSITIVE_NAN = 0x7FC00000, |
| NEGATIVE_NAN = 0xFFC00000; |
| |
| /** |
| * The highest prime number supported by the {@link #nextPrimeNumber(int)} method. |
| * In the current implementation, this value is {@value}. However this limit may |
| * change in any future Apache SIS version. |
| * |
| * <div class="note"><b>Note:</b> |
| * The current value is the highest prime number representable as an unsigned 16 bits integer. |
| * This is enough for current needs because 16 bits prime numbers are sufficient for finding |
| * the divisors of any 32 bits integers.</div> |
| * |
| * @see #nextPrimeNumber(int) |
| */ |
| public static final int HIGHEST_SUPPORTED_PRIME_NUMBER = 65521; |
| |
| /** |
| * Maximal length needed for the {@link #primes} array in order to store prime numbers |
| * from 2 to 32749 (15 bits) or {@value #HIGHEST_SUPPORTED_PRIME_NUMBER} (16 bits). |
| * |
| * @see #primeNumberAt(int) |
| */ |
| static final int PRIMES_LENGTH_15_BITS = 3512, |
| PRIMES_LENGTH_16_BITS = 6542; |
| |
| /** |
| * The sequence of prime numbers computed so far. Will be expanded as needed. |
| * We limit ourself to 16 bits numbers because they are sufficient for computing |
| * divisors of any 32 bits number. |
| * |
| * @see #primeNumberAt(int) |
| */ |
| @SuppressWarnings("VolatileArrayField") // Because we will not modify array content. |
| private static volatile short[] primes = new short[] {2, 3}; |
| |
| /** |
| * Do not allow instantiation of this class. |
| */ |
| private MathFunctions() { |
| } |
| |
| /** |
| * Truncates the given value toward zero. Invoking this method is equivalent to invoking |
| * {@link Math#floor(double)} if the value is positive, or {@link Math#ceil(double)} if |
| * the value is negative. |
| * |
| * @param value the value to truncate. |
| * @return the largest in magnitude (further from zero) integer value which is equals |
| * or less in magnitude than the given value. |
| */ |
| public static double truncate(final double value) { |
| return (doubleToRawLongBits(value) & SIGN_BIT_MASK) == 0 ? Math.floor(value) : Math.ceil(value); |
| } |
| |
| /** |
| * Returns the magnitude of the given vector. This is defined by: |
| * |
| * {@preformat math |
| * sqrt(vector[0]² + vector[1]² + … + vector[length-1]²) |
| * } |
| * |
| * If the given vector contains a NaN value, then the result is NaN. |
| * |
| * <h4>Implementation note</h4> |
| * In the special case where only one element is different than zero, this method |
| * returns directly the {@linkplain Math#abs(double) absolute value} of that element |
| * without computing {@code sqrt(v²)}, in order to avoid rounding error. This special case |
| * has been implemented because this method is often invoked for computing the length of |
| * {@linkplain org.opengis.coverage.grid.RectifiedGrid#getOffsetVectors() offset vectors}, |
| * typically aligned with the axes of a {@linkplain org.opengis.referencing.cs.CartesianCS |
| * Cartesian coordinate system}. |
| * |
| * @param vector the vector for which to compute the magnitude. |
| * @return the magnitude of the given vector as a positive number, of NaN. |
| * |
| * @see Math#hypot(double, double) |
| */ |
| public static double magnitude(final double... vector) { |
| int i = vector.length; |
| |
| // If every elements in the array are zero, returns zero. |
| double v1; |
| do if (i == 0) return 0; |
| while ((v1 = vector[--i]) == 0); |
| |
| // We have found a non-zero element. If it is the only one, returns its absolute value. |
| double v2; |
| do if (i == 0) return abs(v1); |
| while ((v2 = vector[--i]) == 0); |
| |
| // If there is exactly 2 elements, use Math.hypot which is more robust than our algorithm. |
| double v3; |
| do if (i == 0) return Math.hypot(v1, v2); |
| while ((v3 = vector[--i]) == 0); |
| |
| // Usual magnitude computation, but using double-double arithmetic. |
| final DoubleDouble sum = new DoubleDouble(); |
| final DoubleDouble dot = new DoubleDouble(); |
| sum.setToProduct(v1, v1); |
| dot.setToProduct(v2, v2); sum.add(dot); |
| dot.setToProduct(v3, v3); sum.add(dot); |
| while (i != 0) { |
| v1 = vector[--i]; |
| dot.setToProduct(v1, v1); |
| sum.add(dot); |
| } |
| sum.sqrt(); |
| return sum.doubleValue(); |
| } |
| |
| /** |
| * Returns the unbiased exponent used in the representation of a {@code double}, with correction for |
| * sub-normal numbers. This method is related to {@link Math#getExponent(double)} in the following ways: |
| * |
| * <ul> |
| * <li>For NaN and all values equal or greater than {@link Double#MIN_NORMAL} in magnitude (including |
| * infinities), this method returns results that are identical to {@code Math.getExponent(double)}.</li> |
| * <li>For values smaller than {@link Double#MIN_NORMAL} in magnitude (including zero), the correction |
| * for sub-normal numbers results in return values smaller than what {@code Math.getExponent(double)} |
| * would return.</li> |
| * </ul> |
| * |
| * Special cases: |
| * <ul> |
| * <li>If the argument is NaN or infinite, then the result is {@link Double#MAX_EXPONENT} + 1.</li> |
| * <li>If the argument is {@link Double#MAX_VALUE}, then the result is {@value java.lang.Double#MAX_EXPONENT}.</li> |
| * <li>If the argument is {@link Double#MIN_NORMAL}, then the result is {@value java.lang.Double#MIN_EXPONENT}.</li> |
| * <li>If the argument is {@link Double#MIN_VALUE}, then the result is -1074.</li> |
| * <li>If the argument is zero, then the result is -1075.</li> |
| * </ul> |
| * |
| * <h4>Identities</h4> |
| * For any <var>p</var> values in the [-1075 … 1024] range and <var>value</var> = 2<sup>p</sup>: |
| * <ul> |
| * <li><code>getExponent(Math.scalb(1.0, p)) == p</code></li> |
| * <li><code>Math.scalb(1.0, getExponent(value)) == value</code></li> |
| * <li><code>Math.floor({@linkplain #LOG10_2} * getExponent(value)) == Math.floor(Math.log10(value))</code></li> |
| * </ul> |
| * |
| * @param value the value for which to get the exponent. |
| * @return the unbiased exponent, corrected for sub-normal numbers if needed. |
| * Values will be in the [-1075 … 1024] range, inclusive. |
| * |
| * @see Math#getExponent(double) |
| * @see Math#scalb(double, int) |
| * |
| * @since 0.4 |
| */ |
| public static int getExponent(final double value) { |
| final long bits = doubleToRawLongBits(value); |
| int exponent = (int) ((bits >>> SIGNIFICAND_SIZE) & 0x7FFL); |
| if (exponent == 0) { |
| /* |
| * Number is sub-normal: there is no implicit 1 bit before the significand. |
| * We need to search for the position of the first real 1 bit, and fix the |
| * exponent accordingly. Note that numberOfLeadingZeros(…) is relative to |
| * 64 bits while the significand size is only 52 bits. The last term below |
| * is for fixing this difference. |
| */ |
| exponent -= Long.numberOfLeadingZeros(bits & ((1L << SIGNIFICAND_SIZE) - 1)) - (Long.SIZE - SIGNIFICAND_SIZE); |
| } |
| return exponent - Double.MAX_EXPONENT; |
| } |
| |
| /** |
| * Computes the result of {@code base} argument raised to the power given by {@code exponent} argument. |
| * This method computes the same value than {@link Math#pow(double, double)} but using only integer arithmetic. |
| * The result must be representable as a 64 bits integers ({@code long} primitive type), |
| * otherwise an {@link ArithmeticException} is thrown. The result is guaranteed exact, |
| * in contrast to results represented as {@code double} floating point values |
| * which may be approximate for magnitudes greater than 2<sup>52</sup>. |
| * This method may also be faster. |
| * |
| * <div class="note"><b>Implementation note:</b> this method uses |
| * <a href="https://en.wikipedia.org/wiki/Exponentiation_by_squaring">exponentiation by squaring</a> technic.</div> |
| * |
| * The type of the {@code base} argument is {@code long} for convenience, since this method is used in contexts |
| * where relatively large integers are handled. However any value greater than the capacity of {@code int} type |
| * is guaranteed to fail with {@link ArithmeticException} unless {@code exponent} is 0 or 1. |
| * Likewise any {@code exponent} value greater than 62 is guaranteed to fail unless {@code base} is 0 or 1. |
| * |
| * @param base the value to raise to an exponent. |
| * @param exponent the exponent, as zero or positive number. |
| * @return the value <var>base</var><sup><var>exponent</var></sup> as a 64 bits integer. |
| * @throws ArithmeticException if the given exponent is negative, or if the result overflow integer arithmetic. |
| * |
| * @see Math#pow(double, double) |
| * |
| * @since 1.0 |
| */ |
| public static long pow(long base, int exponent) { |
| long result = 1; |
| if (exponent >= 1) { |
| if ((exponent & 1) != 0) { |
| result = base; |
| } |
| while ((exponent >>>= 1) != 0) { |
| base = Math.multiplyExact(base, base); |
| if ((exponent & 1) != 0) { |
| result = Math.multiplyExact(result, base); |
| } |
| } |
| } else if (exponent < 0) { |
| throw new ArithmeticException(Errors.format(Errors.Keys.NegativeArgument_2, "exponent", exponent)); |
| } |
| return result; |
| } |
| |
| /** |
| * Computes 10 raised to the power of <var>x</var>. This method is faster and slightly more accurate |
| * than invoking <code>{@linkplain Math#pow(double, double) Math.pow}(10, x)</code>. |
| * |
| * <div class="note"><b>Note:</b> |
| * This method has been defined because the standard {@code Math.pow(10, x)} method does not always return |
| * the closest IEEE floating point representation. Slight departures (1 or 2 ULP) are often allowed in math |
| * functions for performance reasons. The most accurate calculations are usually not necessary, but the base |
| * 10 is a special case since it is used for scaling axes or formatting human-readable output.</div> |
| * |
| * Special cases: |
| * <ul> |
| * <li>If <var>x</var> is equals or lower than -324, then the result is 0.</li> |
| * <li>If <var>x</var> is equals or greater than 309, then the result is {@linkplain Double#POSITIVE_INFINITY positive infinity}.</li> |
| * <li>If <var>x</var> is in the [0 … 18] range inclusive, then the result is exact.</li> |
| * <li>For all other <var>x</var> values, the result is the closest IEEE 754 approximation.</li> |
| * </ul> |
| * |
| * @param x the exponent. |
| * @return 10 raised to the given exponent. |
| * |
| * @see #pow10(double) |
| * @see #LOG10_2 |
| * @see DecimalFunctions |
| */ |
| public static double pow10(final int x) { |
| return DecimalFunctions.pow10(x); |
| } |
| |
| /** |
| * Computes 10 raised to the power of <var>x</var>. Invoking this method is equivalent to invoking |
| * <code>{@linkplain Math#pow(double, double) Math.pow}(10, x)</code>, but is slightly more accurate |
| * in the special case where the given argument is an integer. |
| * |
| * @param x the exponent. |
| * @return 10 raised to the given exponent. |
| * |
| * @see #pow10(int) |
| * @see Math#pow(double, double) |
| * @see Math#log10(double) |
| */ |
| public static double pow10(final double x) { |
| final int ix = (int) x; |
| if (ix == x) { |
| return DecimalFunctions.pow10(ix); |
| } else { |
| return Math.pow(10, x); |
| } |
| } |
| |
| /** |
| * Returns the inverse hyperbolic sine of the given value. |
| * This is the inverse of the {@link Math#sinh(double)} method. |
| * |
| * @param x the value for which to compute the inverse hyperbolic sine. |
| * @return the inverse hyperbolic sine of the given value. |
| * |
| * @see Math#sinh(double) |
| * |
| * @since 0.6 |
| */ |
| public static double asinh(final double x) { |
| return Math.log(x + sqrt(x*x + 1)); |
| } |
| |
| /** |
| * Returns the inverse hyperbolic cosine of the given value. |
| * This is the inverse of the {@link Math#cosh(double)} method. |
| * |
| * @param x the value for which to compute the inverse hyperbolic cosine. |
| * @return the inverse hyperbolic cosine of the given value. |
| * |
| * @see Math#cosh(double) |
| * |
| * @since 0.6 |
| */ |
| public static double acosh(final double x) { |
| return Math.log(x + sqrt(x*x - 1)); |
| } |
| |
| /** |
| * Returns the inverse hyperbolic tangent of the given value. |
| * This is the inverse of the {@link Math#tanh(double)} method. |
| * The range of input values shall be in the [-1 … 1]. |
| * Special cases: |
| * |
| * <ul> |
| * <li>For <var>x</var> = NaN, this method returns a {@linkplain Double#isNaN(double) NaN} value.</li> |
| * <li>For <var>x</var> = -1, this method returns {@linkplain Double#NEGATIVE_INFINITY negative infinity}.</li> |
| * <li>For <var>x</var> = +1, this method returns {@linkplain Double#POSITIVE_INFINITY positive infinity}.</li> |
| * </ul> |
| * |
| * @param x the value for which to compute the inverse hyperbolic tangent. |
| * @return the inverse hyperbolic tangent of the given value. |
| * |
| * @see Math#tanh(double) |
| */ |
| public static double atanh(final double x) { |
| /* |
| * The classical formulas is log((1+x)/(1-x))/2, but the following is more |
| * accurate if the (1+x)/(1-x) ratio is close to 1, i.e. if x is close to 0. |
| * This is often the case in Apache SIS since x is often a value close to the |
| * Earth excentricity, which is a small value (0 would be a perfect sphere). |
| */ |
| return 0.5 * Math.log1p(2*x / (1-x)); |
| } |
| |
| /** |
| * Returns {@code true} if the given value is positive, <em>excluding</em> negative zero. |
| * Special cases: |
| * |
| * <ul> |
| * <li>If the value is {@code +0.0}, returns {@code true}</li> |
| * <li>If the value is {@code -0.0}, returns <b>{@code false}</b></li> |
| * <li>If the value is {@link Double#isNaN(double) NaN}, returns {@code false}</li> |
| * </ul> |
| * |
| * As seen from the above cases, this method distinguishes positive zero from negative zero. |
| * The handling of zero values is the difference between invoking {@code isPositive(double)} |
| * and testing if (<var>value</var> {@literal >= 0}). |
| * |
| * @param value the value to test. |
| * @return {@code true} if the given value is positive, excluding negative zero. |
| * |
| * @see #isPositiveZero(double) |
| * @see #isNegative(double) |
| */ |
| public static boolean isPositive(final double value) { |
| return (doubleToRawLongBits(value) & SIGN_BIT_MASK) == 0 && !Double.isNaN(value); |
| } |
| |
| /** |
| * Returns {@code true} if the given value is the positive zero ({@code +0.0}). |
| * This method returns {@code false} for the negative zero ({@code -0.0}). |
| * This method is equivalent to the following code, but potentially faster: |
| * |
| * {@preformat java |
| * return (value == 0) && isPositive(value); |
| * } |
| * |
| * @param value the value to test. |
| * @return {@code true} if the given value is +0.0 (not -0.0). |
| * |
| * @see #isPositive(double) |
| * @see #isNegativeZero(double) |
| * |
| * @since 0.4 |
| */ |
| public static boolean isPositiveZero(final double value) { |
| return doubleToRawLongBits(value) == 0L; |
| } |
| |
| /** |
| * Returns {@code true} if the given value is negative, <em>including</em> negative zero. |
| * Special cases: |
| * |
| * <ul> |
| * <li>If the value is {@code +0.0}, returns {@code false}</li> |
| * <li>If the value is {@code -0.0}, returns <b>{@code true}</b></li> |
| * <li>If the value is {@link Double#isNaN(double) NaN}, returns {@code false}</li> |
| * </ul> |
| * |
| * As seen from the above cases, this method distinguishes positive zero from negative zero. |
| * The handling of zero values is the difference between invoking {@code isNegative(double)} |
| * and testing if (<var>value</var> {@literal < 0}). |
| * |
| * @param value the value to test. |
| * @return {@code true} if the given value is negative, including negative zero. |
| * |
| * @see #isNegativeZero(double) |
| * @see #isPositive(double) |
| */ |
| public static boolean isNegative(final double value) { |
| return (doubleToRawLongBits(value) & SIGN_BIT_MASK) != 0 && !Double.isNaN(value); |
| } |
| |
| /** |
| * Returns {@code true} if the given value is the negative zero ({@code -0.0}). |
| * This method returns {@code false} for the positive zero ({@code +0.0}). |
| * This method is equivalent to the following code, but potentially faster: |
| * |
| * {@preformat java |
| * return (value == 0) && isNegative(value); |
| * } |
| * |
| * @param value the value to test. |
| * @return {@code true} if the given value is -0.0 (not +0.0). |
| * |
| * @see #isNegative(double) |
| * @see #isPositiveZero(double) |
| * |
| * @since 0.4 |
| */ |
| public static boolean isNegativeZero(final double value) { |
| return doubleToRawLongBits(value) == SIGN_BIT_MASK; |
| } |
| |
| /** |
| * Returns {@code true} if the given values have the same sign, differentiating positive |
| * and negative zeros. |
| * Special cases: |
| * |
| * <ul> |
| * <li>{@code +0.0} and {@code -0.0} are considered to have opposite sign</li> |
| * <li>If any value is {@link Double#isNaN(double) NaN}, returns {@code false}</li> |
| * </ul> |
| * |
| * @param v1 the first value. |
| * @param v2 the second value, to compare the sign with the first value. |
| * @return {@code true} if the given values are not NaN and have the same sign. |
| * |
| * @see Math#signum(double) |
| */ |
| public static boolean isSameSign(final double v1, final double v2) { |
| return !Double.isNaN(v1) && !Double.isNaN(v2) && |
| ((doubleToRawLongBits(v1) ^ doubleToRawLongBits(v2)) & SIGN_BIT_MASK) == 0; |
| } |
| |
| /** |
| * Returns the first floating-point argument with the sign reversed if the second floating-point |
| * argument is negative. This method is similar to <code>{@linkplain Math#copySign(double,double) |
| * Math.copySign}(value, sign)</code> except that the sign is combined with an <cite>exclusive |
| * or</cite> operation instead than being copied. |
| * |
| * <p>This method makes no guarantee about whether {@code NaN} values are handled as positive |
| * or negative numbers. This is the same policy than {@link Math#copySign(double, double)}.</p> |
| * |
| * @param value the parameter providing the value that may need a sign change. |
| * @param sign the parameter providing the sign to <cite>xor</cite> with the value. |
| * @return the provided value with its sign reversed if the {@code sign} parameter is negative. |
| * |
| * @see Math#copySign(double, double) |
| */ |
| public static double xorSign(final double value, final double sign) { |
| return longBitsToDouble(doubleToRawLongBits(value) ^ |
| (doubleToRawLongBits(sign) & SIGN_BIT_MASK)); |
| } |
| |
| /** |
| * Returns {@code true} if the given values are {@linkplain Float#equals(Object) equal} |
| * or if their difference is not greater than the given threshold. More specifically: |
| * |
| * <ul> |
| * <li>If both values are {@linkplain Float#POSITIVE_INFINITY positive infinity}, or |
| * if both values are {@linkplain Float#NEGATIVE_INFINITY negative infinity}, |
| * then this method returns {@code true}.</li> |
| * <li>If both values {@linkplain Float#isNaN(float) are NaN}, then this method returns {@code true}. |
| * Note that this method does not differentiate the various NaN values.</li> |
| * <li>Otherwise, this method returns the result of the {@code abs(v1 - v2) <= ε} comparison.</li> |
| * </ul> |
| * |
| * @param v1 the first value to compare. |
| * @param v2 the second value to compare. |
| * @param ε the tolerance threshold, which must be positive. |
| * @return {@code true} if both values are equal given the tolerance threshold. |
| */ |
| public static boolean epsilonEqual(final float v1, final float v2, final float ε) { |
| return (abs(v1 - v2) <= ε) || floatToIntBits(v1) == floatToIntBits(v2); |
| } |
| |
| /** |
| * Returns {@code true} if the given values are {@linkplain Double#equals(Object) equal} |
| * or if their difference is not greater than the given threshold. More specifically: |
| * |
| * <ul> |
| * <li>If both values are {@linkplain Double#POSITIVE_INFINITY positive infinity}, or |
| * if both values are {@linkplain Double#NEGATIVE_INFINITY negative infinity}, |
| * then this method returns {@code true}.</li> |
| * <li>If both values {@linkplain Double#isNaN(double) are NaN}, then this method returns {@code true}. |
| * Note that this method does not differentiate the various NaN values.</li> |
| * <li>Otherwise, this method returns the result of the {@code abs(v1 - v2) <= ε} comparison.</li> |
| * </ul> |
| * |
| * @param v1 the first value to compare. |
| * @param v2 the second value to compare. |
| * @param ε the tolerance threshold, which must be positive. |
| * @return {@code true} if both values are equal given the tolerance threshold. |
| */ |
| public static boolean epsilonEqual(final double v1, final double v2, final double ε) { |
| return (abs(v1 - v2) <= ε) || doubleToLongBits(v1) == doubleToLongBits(v2); |
| } |
| |
| /** |
| * Returns a {@linkplain Float#isNaN(float) NaN} number for the specified ordinal value. |
| * Valid NaN numbers in Java can have bit fields in the ranges listed below: |
| * |
| * <ul> |
| * <li>[{@code 0x7F800001} … {@code 0x7FFFFFFF}], with |
| * {@code 0x7FC00000} as the bit fields of the standard {@link Float#NaN} value</li> |
| * <li>[{@code 0xFF800001} … {@code 0xFFFFFFFF}]</li> |
| * </ul> |
| * |
| * Some of those bits, named the <cite>payload</cite>, can be used for storing custom information. |
| * This method maps some of the payload values to each ordinal value. |
| * |
| * <p>This method guarantees that {@code toNanFloat(0)} returns the standard {@link Float#NaN} value. |
| * For all other {@code ordinal} values, the relationship to the payload values is implementation dependent |
| * and may change in any future version of the SIS library. The current implementation restricts the |
| * range of allowed ordinal values to a smaller one than the range of all possible values.</p> |
| * |
| * @param ordinal the NaN ordinal value, from {@value #MIN_NAN_ORDINAL} to {@value #MAX_NAN_ORDINAL} inclusive. |
| * @return one of the legal {@linkplain Float#isNaN(float) NaN} values as a float. |
| * @throws IllegalArgumentException if the specified ordinal is out of range. |
| * |
| * @see Float#intBitsToFloat(int) |
| */ |
| public static float toNanFloat(final int ordinal) throws IllegalArgumentException { |
| ArgumentChecks.ensureBetween("ordinal", MIN_NAN_ORDINAL, MAX_NAN_ORDINAL, ordinal); |
| int bits = (ordinal >= 0) ? ordinal : ~ordinal; |
| bits = (bits + POSITIVE_NAN) | (ordinal & Integer.MIN_VALUE); |
| assert Integer.compareUnsigned(bits, ordinal >= 0 ? POSITIVE_NAN : NEGATIVE_NAN) >= 0 : ordinal; |
| final float value = intBitsToFloat(bits); |
| assert Float.isNaN(value) && toNanOrdinal(value) == ordinal : ordinal; |
| return value; |
| } |
| |
| /** |
| * Returns the ordinal value of the given NaN number. |
| * This method is the converse of {@link #toNanFloat(int)}. |
| * |
| * <p>If the given float is the standard {@link Float#NaN} value, then this method returns 0. |
| * For all other values, the relationship between the float payload and the returned ordinal |
| * is implementation dependent and may change in any future Apache SIS version.</p> |
| * |
| * @param value the value from which to get the NaN ordinal value. |
| * @return the NaN ordinal value of the given floating point value. |
| * @throws IllegalArgumentException if the given value is not a NaN value, |
| * or does not use a supported bits pattern. |
| */ |
| public static int toNanOrdinal(final float value) throws IllegalArgumentException { |
| final int bits = floatToRawIntBits(value); |
| int ordinal = (bits & Integer.MAX_VALUE) - POSITIVE_NAN; |
| if (bits < 0) ordinal = ~ordinal; |
| if (ordinal >= MIN_NAN_ORDINAL && ordinal <= MAX_NAN_ORDINAL) { |
| return ordinal; |
| } |
| final short resourceKey; |
| final Object obj; |
| if (Float.isNaN(value)) { |
| resourceKey = Errors.Keys.IllegalBitsPattern_1; |
| obj = Integer.toHexString(bits); |
| } else { |
| resourceKey = Errors.Keys.IllegalArgumentValue_2; |
| obj = new Object[] {"value", value}; |
| } |
| throw new IllegalArgumentException(Errors.format(resourceKey, obj)); |
| } |
| |
| /** |
| * Converts two long bits values containing a IEEE 754 quadruple precision floating point number |
| * to a double precision floating point number. About 17 decimal digits of precision may be lost |
| * due to the {@code double} type having only half the capacity of quadruple precision type. |
| * |
| * <p>Some quadruple precision values can not be represented in double precision and are mapped |
| * to {@code double} values as below:</p> |
| * <ul> |
| * <li>Values having a magnitude less than {@link Double#MIN_VALUE} are mapped to |
| * positive or negative zero.</li> |
| * <li>Values having a magnitude greater than {@link Double#MAX_VALUE} are mapped to |
| * {@link Double#POSITIVE_INFINITY} or {@link Double#NEGATIVE_INFINITY}.</li> |
| * <li>All NaN values are currently collapsed to the single "canonical" {@link Double#NaN} value |
| * (this policy may be revisited in future SIS version).</li> |
| * </ul> |
| * |
| * @param l0 upper part of the quadruple precision floating point number. |
| * @param l1 lower part of the quadruple precision floating point number. |
| * @return double precision approximation. |
| * |
| * @see <a href="https://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format">Quadruple-precision floating-point format on Wikipedia</a> |
| * |
| * @since 0.7 |
| */ |
| public static double quadrupleToDouble(long l0, long l1) { |
| // Build double |
| long sig = (l0 & 0x8000000000000000L); |
| long exp = (l0 & 0x7FFF000000000000L) >> 48; |
| l0 = (l0 & 0x0000FFFFFFFFFFFFL); |
| if (exp == 0) { |
| /* |
| * Subnormal number. |
| * Since we convert them to double precision, subnormal numbers can not be represented |
| * as they are smaller than Double.MIN_VALUE. We map them to zero preserving the sign. |
| */ |
| return Double.longBitsToDouble(sig); |
| } |
| if (exp == 0x7FFF) { |
| /* |
| * NaN of infinite number. |
| * Mantissa with all bits at 0 is used for infinite. |
| * This is the only special number that we can preserve. |
| */ |
| if (l0 == 0 && l1 == 0) { |
| return Double.longBitsToDouble(sig | 0x7FF0000000000000L); |
| } |
| /* |
| * Other NaN values might have a meaning (e.g. NaN(1) = forest, NaN(2) = lake, etc.) |
| * See above toNanFloat(int) and toNaNOrdinal(float) methods. When truncating the value we |
| * might change the meaning, which could cause several issues later. Therefor we conservatively |
| * collapse all NaNs to the default NaN for now (this may be revisited in a future SIS version). |
| */ |
| return Double.NaN; |
| } |
| exp -= (16383 - 1023); //change from 15 bias to 11 bias |
| // Check cases where mantissa excess what double can support |
| if (exp < 0) return Double.NEGATIVE_INFINITY; |
| if (exp > 2046) return Double.POSITIVE_INFINITY; |
| |
| return Double.longBitsToDouble(sig | (exp << 52) | (l0 << 4) | (l1 >>> 60)); |
| } |
| |
| /** |
| * Returns the <var>i</var><sup>th</sup> prime number. |
| * This method returns (2, 3, 5, 7, 11, …) for index (0, 1, 2, 3, 4, …). |
| * |
| * @param index the prime number index, starting at index 0 for prime number 2. |
| * @return the prime number at the specified index. |
| * @throws IndexOutOfBoundsException if the specified index is too large. |
| * |
| * @see java.math.BigInteger#isProbablePrime(int) |
| */ |
| static int primeNumberAt(final int index) throws IndexOutOfBoundsException { |
| ArgumentChecks.ensureValidIndex(PRIMES_LENGTH_16_BITS, index); |
| short[] primes = MathFunctions.primes; |
| if (index >= primes.length) { |
| synchronized (MathFunctions.class) { |
| primes = MathFunctions.primes; |
| if (index >= primes.length) { |
| int i = primes.length; |
| int n = Short.toUnsignedInt(primes[i - 1]); |
| // Compute by block of 16 values, for reducing the amount of array resize. |
| primes = Arrays.copyOf(primes, Math.min((index | 0xF) + 1, PRIMES_LENGTH_16_BITS)); |
| do { |
| testNextNumber: while (true) { // Simulate a "goto" statement (usually not recommanded...) |
| final int stopAt = (int) sqrt(n += 2); |
| int prime; |
| int j = 0; |
| do { |
| prime = Short.toUnsignedInt(primes[++j]); |
| if (n % prime == 0) { |
| continue testNextNumber; |
| } |
| } while (prime <= stopAt); |
| primes[i] = (short) n; |
| break; |
| } |
| } while (++i < primes.length); |
| MathFunctions.primes = primes; |
| } |
| } |
| } |
| return Short.toUnsignedInt(primes[index]); |
| } |
| |
| /** |
| * Returns the first prime number equals or greater than the given value. |
| * Current implementation accepts only values in the |
| * [2 … {@value #HIGHEST_SUPPORTED_PRIME_NUMBER}] range. |
| * |
| * @param number the number for which to find the next prime. |
| * @return the given number if it is a prime number, or the next prime number otherwise. |
| * @throws IllegalArgumentException if the given value is outside the supported range. |
| * |
| * @see java.math.BigInteger#isProbablePrime(int) |
| */ |
| public static int nextPrimeNumber(final int number) throws IllegalArgumentException { |
| ArgumentChecks.ensureBetween("number", 2, HIGHEST_SUPPORTED_PRIME_NUMBER, number); |
| final short[] primes = MathFunctions.primes; |
| int lower = 0; |
| int upper = Math.min(PRIMES_LENGTH_15_BITS, primes.length); |
| if (number > Short.MAX_VALUE) { |
| lower = upper; |
| upper = primes.length; |
| } |
| int i = Arrays.binarySearch(primes, lower, upper, (short) number); |
| if (i < 0) { |
| i = ~i; |
| if (i >= primes.length) { |
| int p; |
| do p = primeNumberAt(i++); |
| while (p < number); |
| return p; |
| } |
| } |
| return Short.toUnsignedInt(primes[i]); |
| } |
| |
| /** |
| * Returns the divisors of the specified number as positive integers. For any value other |
| * than {@code O} (which returns an empty array), the first element in the returned array |
| * is always {@code 1} and the last element is always the absolute value of {@code number}. |
| * |
| * @param number the number for which to compute the divisors. |
| * @return the divisors in strictly increasing order. |
| */ |
| public static int[] divisors(int number) { |
| if (number == 0) { |
| return ArraysExt.EMPTY_INT; |
| } |
| number = abs(number); |
| int[] divisors = new int[16]; |
| divisors[0] = 1; |
| int count = 1; |
| /* |
| * Searches for the first divisors among the prime numbers. We stop the search at the |
| * square root of 'n' because every values above that point can be inferred from the |
| * values before that point, i.e. if n=p1*p2 and p2 is greater than 'sqrt', than p1 |
| * most be lower than 'sqrt'. |
| */ |
| final int sqrt = (int) sqrt(number); // Really want rounding toward 0. |
| for (int p,i=0; (p=primeNumberAt(i)) <= sqrt; i++) { |
| if (number % p == 0) { |
| if (count == divisors.length) { |
| divisors = Arrays.copyOf(divisors, count*2); |
| } |
| divisors[count++] = p; |
| } |
| } |
| /* |
| * Completes the divisors past 'sqrt'. The numbers added here may or may not be prime |
| * numbers. Side note: checking that they are prime numbers would be costly, but this |
| * algorithm doesn't need that. |
| */ |
| int source = count; |
| if (count*2 > divisors.length) { |
| divisors = Arrays.copyOf(divisors, count*2); |
| } |
| int d1 = divisors[--source]; |
| int d2 = number / d1; |
| if (d1 != d2) { |
| divisors[count++] = d2; |
| } |
| while (--source >= 0) { |
| divisors[count++] = number / divisors[source]; |
| } |
| /* |
| * Checks the products of divisors found so far. For example if 2 and 3 are divisors, |
| * checks if 6 is a divisor as well. The products found will themself be used for |
| * computing new products. |
| */ |
| for (int i=1; i<count; i++) { |
| d1 = divisors[i]; |
| for (int j=i; j<count; j++) { |
| d2 = d1 * divisors[j]; |
| if (number % d2 == 0) { |
| int p = Arrays.binarySearch(divisors, j, count, d2); |
| if (p < 0) { |
| p = ~p; // tild (~) operator, not minus |
| if (count == divisors.length) { |
| divisors = Arrays.copyOf(divisors, count*2); |
| } |
| System.arraycopy(divisors, p, divisors, p+1, count-p); |
| divisors[p] = d2; |
| count++; |
| } |
| } |
| } |
| } |
| divisors = ArraysExt.resize(divisors, count); |
| assert ArraysExt.isSorted(divisors, true); |
| return divisors; |
| } |
| |
| /** |
| * Returns the positive divisors which are common to all the specified numbers. |
| * The returned array always starts with value 1, unless the given value is 0 |
| * in which case this method returns an empty array. |
| * |
| * @param numbers the numbers for which to compute the divisors, in any order. |
| * @return the divisors common to all the given numbers, in strictly increasing order. |
| */ |
| public static int[] commonDivisors(final int... numbers) { |
| if (numbers.length == 0) { |
| return ArraysExt.EMPTY_INT; |
| } |
| /* |
| * Get the smallest value. We will compute the divisors only for this value, |
| * since we know that any value greater that the minimal value can not be a |
| * common divisor. |
| */ |
| int minValue = Integer.MAX_VALUE; |
| for (int i=0; i<numbers.length; i++) { |
| final int n = abs(numbers[i]); |
| if (n <= minValue) { |
| minValue = n; |
| } |
| } |
| int[] divisors = divisors(minValue); |
| /* |
| * Tests if the divisors we just found are also divisors of all other numbers. |
| * Removes those which are not. |
| */ |
| int count = divisors.length; |
| for (int i=0; i<numbers.length; i++) { |
| final int n = abs(numbers[i]); |
| if (n != minValue) { |
| for (int j=count; --j>0;) { // Do not test j==0, since divisors[0] == 1. |
| if (n % divisors[j] != 0) { |
| System.arraycopy(divisors, j+1, divisors, j, --count - j); |
| } |
| } |
| } |
| } |
| return ArraysExt.resize(divisors, count); |
| } |
| |
| /** |
| * Returns the real (non-complex) roots of a polynomial equation having the given coefficients. |
| * This method returns the <var>x</var> values for which <var>y</var>=0 in the following equation: |
| * |
| * <blockquote><var>y</var> = |
| * <var>c<sub>0</sub></var> + |
| * <var>c<sub>1</sub></var>⋅<var>x</var> + |
| * <var>c<sub>2</sub></var>⋅<var>x</var><sup>2</sup> + |
| * <var>c<sub>3</sub></var>⋅<var>x</var><sup>3</sup> + … + |
| * <var>c<sub>n</sub></var>⋅<var>x</var><sup>n</sup> |
| * </blockquote> |
| * |
| * Current implementation can resolve polynomials described by a maximum of 5 coefficients, ignoring |
| * leading and trailing zeros. They correspond to linear, quadratic, cubic and quartic polynomials. |
| * |
| * @param coefficients the <var>c<sub>0</sub></var>, <var>c<sub>1</sub></var>, <var>c<sub>2</sub></var>, … |
| * <var>c<sub>n</sub></var> coefficients, in that order. |
| * @return the non-complex roots, or an empty array if none. |
| * @throws UnsupportedOperationException if given arguments contain more non-zero coefficients than this method can handle. |
| * |
| * @see java.awt.geom.QuadCurve2D#solveQuadratic(double[]) |
| * @see java.awt.geom.CubicCurve2D#solveCubic(double[]) |
| * |
| * @since 1.0 |
| */ |
| public static double[] polynomialRoots(final double... coefficients) { |
| int upper = coefficients.length; |
| while (upper > 0) { |
| double a = coefficients[--upper]; // Coefficient of the term with highest exponent. |
| if (a == 0) { |
| continue; // Search the last non-zero coefficient. |
| } |
| double c; // Coefficient of the term with lowest exponent. |
| int lower = 0; |
| while ((c = coefficients[lower]) == 0) lower++; |
| switch (upper - lower) { |
| /* |
| * c = 0 |
| * |
| * Can not compute x. We could return an arbitrary value if c = 0, but we rather return |
| * an empty array for keeping the number of root equals to the highest exponent. |
| */ |
| case 0: { |
| break; |
| } |
| /* |
| * ax + c = 0 → x = -c/a |
| */ |
| case 1: { |
| final double x = -c / a; |
| if (Double.isNaN(x)) break; |
| return new double[] {x}; |
| } |
| /* |
| * ax² + bx + c = 0 → x = (-b ± √[b² - 4ac]) / (2a) |
| * |
| * Above equation is numerically unstable. More stable algorithm is given |
| * by Numerical Recipes §5.6 (quadratic equation), which is applied below. |
| */ |
| case 2: { |
| final double b = coefficients[lower + 1]; |
| final double q = -0.5 * (b + Math.copySign(sqrt(b*b - 4*a*c), b)); |
| final double x1 = q/a; |
| final double x2 = c/q; |
| if (Double.isNaN(x1) && Double.isNaN(x2)) break; |
| return (x1 != x2) ? new double[] {x1, x2} : new double[] {x1}; |
| } |
| /* |
| * x³ + ax² + bx + c = 0 |
| * |
| * Numerical Recipes §5.6 (cubic equation) is applied below. |
| * Solution usually have either 1 or 3 roots. |
| */ |
| case 3: { |
| return refineRoots(coefficients, solveCubic( |
| coefficients[lower + 2] / a, |
| coefficients[lower + 1] / a, |
| c / a, false)); |
| } |
| /* |
| * x⁴ + ax³ + bx² + cx + d = 0 |
| * |
| * https://dlmf.nist.gov/1.11 in "Quartic equations" section. |
| * This algorithm reduces the equation to a cubic equation. |
| */ |
| case 4: { |
| double b,d; |
| d = c / a; |
| c = coefficients[lower + 1] / a; |
| b = coefficients[lower + 2] / a; |
| a = coefficients[lower + 3] / a; |
| final double a2 = a*a; |
| final double p = -3./8 * (a2) + b; |
| final double q = 1./8 * (a2*a) - 1./2 * (a*b) + c; |
| final double r = -3./256 * (a2*a2) + 1./16 * (a2*b) - 1./4 * (a*c) + d; |
| final double[] roots = solveCubic(-2*p, p*p-4*r, q*q, true); |
| if (roots.length != 4) break; |
| for (int i=0; i<3; i++) { |
| roots[i] = sqrt(-roots[i]); |
| } |
| if (isPositive(q)) { |
| for (int i=0; i<3; i++) { |
| roots[i] = -roots[i]; |
| } |
| } |
| final double α = roots[0], β = roots[1], γ = roots[2]; |
| final double s = α + β + γ; |
| if (Double.isNaN(s)) break; |
| roots[0] = s/2 - (a /= 4); |
| roots[1] = (+α - β - γ)/2 - a; |
| roots[2] = (-α + β - γ)/2 - a; |
| roots[3] = (-α - β + γ)/2 - a; |
| return refineRoots(coefficients, removeDuplicated(roots)); |
| } |
| default: { |
| throw new UnsupportedOperationException(); |
| } |
| } |
| break; |
| } |
| return ArraysExt.EMPTY_DOUBLE; |
| } |
| |
| /** |
| * Solves cubic equation x³ + ax² + bx + c = 0. The solution before simplification has either 1 or 3 roots. |
| * The {@code quartic} argument specifies whether this cubic equation is used as a step for solving a quartic equation: |
| * |
| * <ul> |
| * <li>If {@code true}, then we are interested only in the 3 roots solution and we do not check for duplicated values. |
| * The length of returned array is 4 for allowing the caller to reuse the same array.</li> |
| * <li>If {@code false}, then this method may simplify the 3 roots to 2 roots if two of them are equal, |
| * or may return the 1 root solution. The length of returned array is the number of roots.</li> |
| * </ul> |
| */ |
| private static double[] solveCubic(double a, double b, double c, final boolean quartic) { |
| final double Q = (a*a - 3*b) / 9; // Q from Numerical Recipes 5.6.10. |
| final double R = (a*(a*a - 4.5*b) + 13.5*c) / 27; // R from Numerical Recipes 5.6.10. |
| final double Q3 = Q*Q*Q; |
| final double R2 = R*R; |
| a /= 3; // Last term of Numerical Recipes 5.6.12, 17 and 18. |
| if (R2 < Q3) { |
| /* |
| * Numerical Recipes 5.6.11 and 5.6.12 uses acos(R/sqrt(Q³)). It is possible to rewrite as |
| * atan2(sqrt(Q3 - R2), R) using the cos(θ) = 1/√[1 + tan²θ] trigonometric identity, but |
| * this substitution seems to decrease accuracy instead of increasing it in our tests. |
| */ |
| b = Math.acos(R/sqrt(Q3)) / 3; // θ from Numerical recipes 5.6.11, then b = θ/3. |
| c = -2 * sqrt(Q); // First part of Numerical Recipes 5.6.12. |
| double[] roots = new double[quartic ? 4 : 3]; |
| roots[2] = c*Math.cos(b - 2*Math.PI/3) - a; // TODO: try Math.fma with JDK9. |
| roots[1] = c*Math.cos(b + 2*Math.PI/3) - a; |
| roots[0] = c*Math.cos(b) - a; |
| if (!quartic) { |
| roots = removeDuplicated(roots); |
| } |
| return roots; |
| } |
| if (!quartic) { |
| b = -Math.copySign(Math.cbrt(abs(R) + sqrt(R2 - Q3)), R); // A from Numerical Recipes 5.6.15. |
| final double x = (b == 0 ? 0 : b + Q/b) - a; |
| if (!Double.isNaN(x)) { |
| return new double[] {x}; |
| } |
| } |
| return ArraysExt.EMPTY_DOUBLE; |
| } |
| |
| /** |
| * Remove duplicated values in the given array. This method is okay only for very small arrays (3 or 4 elements). |
| * Duplicated values should be very rare and occur mostly as a consequence of rounding errors while computing the |
| * roots of polynomial equations. Because if the algebraic solution has less roots than what we would expect from |
| * the largest exponent (for example ax² + bx = 0 has only one root instead of two), then {@link #polynomialRoots} |
| * should have reduced the equation to a lower degrees (ax + b = 0 in above example), in which case there is no |
| * duplicated roots to remove. |
| */ |
| private static double[] removeDuplicated(double[] roots) { |
| int i = 1; |
| next: while (i < roots.length) { |
| for (int j=i; --j >= 0;) { |
| if (roots[j] == roots[i]) { |
| roots = ArraysExt.remove(roots, i, 1); |
| continue next; |
| } |
| } |
| i++; |
| } |
| return roots; |
| } |
| |
| /** |
| * Tries to improves accuracy of polynomial roots by applying small displacements |
| * to the <var>x</var> values using ∂y/∂x derivative around those values. |
| * |
| * <div class="note"><b>Purpose:</b> |
| * this refinement is significant in a {@link org.apache.sis.referencing.GeodesicsOnEllipsoid} |
| * test checking the value of an μ(x²,y²) function.</div> |
| * |
| * @param coefficients the user-specified coefficients. |
| * @param roots the roots. This array will be modified in place. |
| * @return {@code roots}. |
| */ |
| private static double[] refineRoots(final double[] coefficients, final double[] roots) { |
| for (int i=0; i < roots.length; i++) { |
| double ymin = Double.POSITIVE_INFINITY; |
| double x = roots[i]; |
| double dx; |
| do { |
| double px = 1; // Power of x: 1, x¹, x², x³, … |
| double dy = 0; // First derivative of polynomial at x. |
| double y = coefficients[0]; // Value of polynomial at x. |
| double ey = 0, edy = 0; // Error terms for Kahan summation algorithm. |
| for (int j=1; j<coefficients.length; j++) { |
| final double c = coefficients[j]; |
| double s; |
| s = c * (px * i) + edy; edy = s + (dy - (dy += s)); // Kahan summation of dy. |
| s = c * (px *= x) + ey; ey = s + ( y - ( y += s)); // Kahan summation of y. |
| } |
| if (!(ymin > (ymin = abs(y)))) break; // If result not better than previous result, stop. |
| roots[i] = x; |
| dx = y/dy; |
| } while (x != (x -= dx) && Double.isFinite(x)); |
| } |
| return roots; |
| } |
| } |
