Google Git
Sign in
apache / sis / 76bddb705ed107711db475c6f348b79748ab850d / . / endorsed / src / org.apache.sis.util / main / org / apache / sis / util / privy / DoubleDouble.java
blob: 15f522bcde882fb45a695577373cf716bcd8f2a4 [file] [log] [blame]
/*
* 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.util.privy;
import java.util.Arrays;
import java.time.Duration;
import java.math.BigInteger;
import java.math.BigDecimal;
import java.math.MathContext;
import org.apache.sis.math.Fraction;
import org.apache.sis.math.MathFunctions;
import org.apache.sis.math.DecimalFunctions;
import org.apache.sis.system.Configuration;
import org.apache.sis.util.Debug;
import static org.apache.sis.pending.jdk.JDK19.DOUBLE_PRECISION;
/**
* Basic arithmetic methods for extended precision numbers using the <i>double-double</i> algorithm.
* This class implements some of the methods published in the following paper:
*
* <ul>
* <li>Yozo Hida, Xiaoye S. Li, David H. Bailey.
* <a href="http://web.mit.edu/tabbott/Public/quaddouble-debian/qd-2.3.4-old/docs/qd.pdf">Library
* for Double-Double and Quad-Double arithmetic</a>, 2007.</li>
* <li>Jonathan R. Shewchuk. Adaptive precision floating-point arithmetic and fast robust geometric predicates.
* Discrete &amp; Computational Geometry, 18(3):305–363, 1997.</li>
* </ul>
*
* {@code DoubleDouble} is used as an alternative to {@link java.math.BigDecimal} when we do not need arbitrary
* precision, we do not want to convert from base 2 to base 10, we need support for NaN and infinities, we want
* more compact storage and better performance. {@code DoubleDouble} can be converted to {@code BigDecimal} as
* below:
*
* {@snippet lang="java" :
* BigDecimal decimal = new BigDecimal(dd.value).add(new BigDecimal(dd.error));
* }
*
* {@code DoubleDouble} is a <em>value object</em>, immutable and thread-safe.
*
* @author Martin Desruisseaux (Geomatys)
*
* @see <a href="https://en.wikipedia.org/wiki/Double-double_%28arithmetic%29#Double-double_arithmetic">Wikipedia: Double-double arithmetic</a>
*/
public final class DoubleDouble extends Number implements Comparable<DoubleDouble> {
/**
* For cross-version compatibility.
*/
private static final long serialVersionUID = -7602414219228638550L;
/**
* {@code true} for disabling the extended precision. This variable should always be {@code false},
* except for testing purpose. If set to {@code true}, then all double-double arithmetic operations
* are immediately followed by a clearing of {@link DoubleDouble#error}. The result should then be
* identical to computation performed using the normal {@code double} arithmetic.
*
* <p>Because this flag is static final, all expressions of the form {@code if (DISABLED)} should be
* omitted by the compiler from the class files in normal operations.</p>
*
* <p>Setting this flag to {@code true} causes some JUnit tests to fail, which is expected. The main
* purpose of this flag is to allow {@code org.apache.sis.referencing.operation.matrix.MatrixTestCase}
* to perform strict comparisons of matrix operation results with JAMA, which is taken as the reference
* implementation. Since JAMA uses {@code double} arithmetic, SIS needs to disable double-double
* arithmetic if the results are to be compared for strict equality.</p>
*/
@Debug
public static final boolean DISABLED = false;
/**
* A margin in number of bits for determining if an error term should be considered as zero.
* When computing <var>a</var> - <var>b</var> as a double-double (106 significand bits) value,
* if the number of non-zero significand bits is equal or lower than {@code ZERO_THRESHOLD+1},
* consider the result as zero.
*/
@Configuration
private static final int ZERO_THRESHOLD = 2;
/**
* Predefined constants frequently used in SIS, sorted in increasing order. This table contains only
* constants that cannot be inferred by {@link DecimalFunctions#deltaForDoubleToDecimal(double)},
* for example some transcendental values.
*
* <p>Elements in this array shall be sorted in strictly increasing order.
* For any value at index <var>i</var>, the associated error is {@code ERRORS[i]}.
*
* @see #errorForWellKnownValue(double)
*/
private static final double[] VALUES = {
/*
* Some of the following constants have more fraction digits than necessary. We declare the extra
* digits for documentation purpose, and in order to have identical content than DoubleDoubleTest
* so that a plain copy-and-paste can be performed between those two classes.
*/
0.000004848136811095359935899141023579480, // Arc-second to radians
0.0002777777777777777777777777777777778, // Second to degrees
0.002777777777777777777777777777777778, // 1/360°
0.01666666666666666666666666666666667, // Minute to degrees
0.01745329251994329576923690768488613, // Degree to radians
0.785398163397448309615660845819876, // π/4
1.111111111111111111111111111111111, // Grad to degrees
1.414213562373095048801688724209698, // √2
1.570796326794896619231321691639751, // π/2
2.356194490192344928846982537459627, // π * 3/4
3.14159265358979323846264338327950, // π
6.28318530717958647692528676655901, // 2π
57.2957795130823208767981548141052 // Radians to degrees
};
/**
* The errors associated to the values in the {@link #VALUES} array.
*
* <p>Tips:</p>
* <ul>
* <li>To compute a new value in this array, just put zero and execute
* {@code DoubleDoubleTest.testErrorForWellKnownValue()}.
* The error message will give the expected value.</li>
* <li>If a computed value is zero, then there is no point to create an entry
* in the {@code (VALUES, ERRORS)} arrays for that value.</li>
* </ul>
*/
private static final double[] ERRORS = {
/* 0.000004… */ 9.320078015422868E-23,
/* 0.000266… */ 2.4093381610788987E-22,
/* 0.002666… */ -1.0601087908747154E-19,
/* 0.016666… */ 2.312964634635743E-19,
/* 0.017453… */ 2.9486522708701687E-19,
/* 0.785398… */ 3.061616997868383E-17,
/* 1.111111… */ -4.9343245538895844E-17,
/* 1.414213… */ -9.667293313452913E-17,
/* 1.570796… */ 6.123233995736766E-17,
/* 2.356194… */ 9.184850993605148E-17,
/* 3.141592… */ 1.2246467991473532E-16,
/* 6.283185… */ 2.4492935982947064E-16,
/* 57.295779… */ -1.9878495670576283E-15
};
/**
* A {@code DoubleDouble} instance for the 0 value.
*
* @see #isZero()
*/
public static final DoubleDouble ZERO = new DoubleDouble(0, 0);
/**
* A {@code DoubleDouble} instance for the 1 value.
*/
public static final DoubleDouble ONE = new DoubleDouble(1, 0);
/**
* A {@code DoubleDouble} instance for the π value.
*/
public static final DoubleDouble PI =
new DoubleDouble(3.14159265358979323846264338327950, 1.2246467991473532E-16);
/**
* A {@code DoubleDouble} instance for the conversion factor from radians to angular degrees.
*/
public static final DoubleDouble RADIANS_TO_DEGREES =
new DoubleDouble(57.2957795130823208767981548141052, -1.9878495670576283E-15);
/**
* A {@code DoubleDouble} instance for the conversion factor from angular degrees to radians.
*/
public static final DoubleDouble DEGREES_TO_RADIANS =
new DoubleDouble(0.01745329251994329576923690768488613, 2.9486522708701687E-19);
/**
* A {@code DoubleDouble} instance for the conversion factor from arc-seconds to radians.
*/
public static final DoubleDouble SECONDS_TO_RADIANS =
new DoubleDouble(0.000004848136811095359935899141023579480, 9.320078015422868E-23);
/**
* A {@code DoubleDouble} instance for the NaN value.
*
* @see #isNaN()
*/
public static final DoubleDouble NaN = new DoubleDouble(Double.NaN, Double.NaN);
/**
* The main value, minus the {@link #error}.
*/
public final double value;
/**
* The error that shall be added to the main {@link #value} in order to get the
* <q>real</q> (actually <q>the most accurate that we can</q>) value.
*/
public final double error;
/**
* Creates a new instance for the given value and error.
* It is caller's responsibility to ensure that the (value, error) pair is normalized.
*
* @param value the value.
* @param error the error.
*/
private DoubleDouble(final double value, final double error) {
this.value = value;
this.error = DISABLED ? 0 : error;
assert !(Math.abs(error) >= Math.ulp(value)) : this; // Use ! for being tolerant to NaN.
}
/**
* Returns an instance for the given number. If the given number is an instance of
* {@code DoubleDouble}, {@link BigDecimal}, {@link BigInteger} or {@link Fraction},
* then the error term will be taken in account.
*
* @param value the value, or {@code null}.
* @param decimal whether {@code float} and {@code double} values were intended to be exact in base 10.
* @return the value as a double-double number, or {@code null} if the given number was null.
*/
public static DoubleDouble of(Number value, final boolean decimal) {
if (value == null) {
return null;
}
if (value instanceof DoubleDouble) {
return (DoubleDouble) value;
}
if (value instanceof Fraction) {
final Fraction f = (Fraction) value;
return new DoubleDouble(f.numerator, 0)
.divide(f.denominator);
}
final double v, error;
if (value instanceof Float) {
final float f = (Float) value;
value = decimal ? DecimalFunctions.floatToDouble(f) : f;
} else if (value instanceof BigInteger) {
value = new BigDecimal((BigInteger) value, MathContext.DECIMAL128);
}
v = value.doubleValue();
if (value instanceof Integer) {
error = 0;
} else if (value instanceof Long) {
error = value.longValue() - (long) v; // Need rounding toward zero.
} else if (value instanceof BigDecimal) {
// Really need new BigDecimal(value) below, not BigDecimal.valueOf(value).
error = ((BigDecimal) value).subtract(new BigDecimal(v), MathContext.DECIMAL64).doubleValue();
} else {
error = decimal ? errorForWellKnownValue(v) : 0;
}
return new DoubleDouble(v, error);
}
/**
* Returns an instance for the given duration in nanoseconds.
*
* @param value the duration to convert.
* @return the given duration, in nanoseconds.
*/
public static DoubleDouble of(final Duration value) {
return of(value.getSeconds()).multiply(Constants.NANOS_PER_SECOND).add(value.getNano());
}
/**
* Returns an instance for the given integer.
*
* @param value the integer value to wrap in a {@code DoubleDouble}.
* @return the value as a double-double number.
*/
public static DoubleDouble of(final int value) {
return new DoubleDouble(value, 0);
}
/**
* Returns an instance for the given long integer.
*
* @param value the long integer value to wrap in a {@code DoubleDouble}.
* @return the value as a double-double number.
*/
public static DoubleDouble of(final long value) {
final double f = value;
return new DoubleDouble(f, value - (long) f); // Rounding toward zero.
}
/**
* Returns an instance for the given value.
* If {@code decimal} is {@code true}, then an error term is inferred for well-known values.
* {@code decimal} should be {@code false} when the value has been computed using transcendental functions
* (cosine, logarithm, <i>etc.</i>), in which case there is no way we can infer a meaningful error term.
* Should also be {@code false} (for performance reason) when the value is an exact representation in base 2.
*
* @param value the value.
* @param decimal whether the value was intended to be exact in base 10.
* @return the value as a double-double number.
*/
public static DoubleDouble of(final double value, final boolean decimal) {
return new DoubleDouble(value, decimal ? errorForWellKnownValue(value) : 0);
}
/**
* Returns an instance for the given value and error.
* It is caller's responsibility to ensure the {@code error} term is less than 1 ULP of {@code value}.
*
* @param value the value.
* @param error the error.
* @return the value as a double-double number.
*/
public static DoubleDouble of(final double value, final double error) {
return new DoubleDouble(value, error);
}
/**
* Return value + error. The result should be identical to {@link #value},
* but we nevertheless do the sum as a safety.
*
* @return {@link #value} + {@link #error}.
*/
@Override public double doubleValue() {return value + error;}
@Override public float floatValue() {return (float) doubleValue();}
@Override public int intValue() {return Numerics.clamp(longValue());}
@Override public long longValue() {return Numerics.saturatingAdd(Math.round(value), (long) error);}
/*
* Do not override `shortValue()` and `byteValue()` in order to keep a behavior
* consistent with all `Number` subclasses provided in the standard JDK: first
* a narrowing conversion to `int` followed by discarding the high order bits.
* Note than even a direct `(short) value` cast implicitly does above steps.
*/
/**
* Suggests an {@link #error} for the given value. The {@code DoubleDouble} class contains a hard-coded list
* of some frequently used constants, for example for various factors of π. If the given value matches exactly
* one of those constants, then its error term is returned. Otherwise this method assumes that the given value
* is defined in base 10 (e.g. many unit conversion factors) and tries to compute an error term with
* {@link DecimalFunctions#deltaForDoubleToDecimal(double)}.
*
* <h4>Rational</h4>
* SIS often creates matrices for unit conversions, and most conversion factors are defined precisely in base 10.
* For example, the conversion from feet to metres is defined by a factor of exactly 0.3048, which cannot be
* represented precisely as a {@code double}. Consequently, if a value of 0.3048 is given, we can assume that
* the intent was to provide the "feet to metres" conversion factor and complete the double-double instance
* accordingly.
*
* @param value the value for which to get this error.
* @return the error for the given value, or 0 if unknown. In the latter case,
* the base 2 representation of the given value is assumed to be accurate enough.
*/
static double errorForWellKnownValue(final double value) {
if (DISABLED) return 0;
final int i = Arrays.binarySearch(VALUES, Math.abs(value));
final double error;
if (i >= 0) {
error = MathFunctions.xorSign(ERRORS[i], value);
} else {
final double delta = DecimalFunctions.deltaForDoubleToDecimal(value);
error = Double.isNaN(delta) ? 0 : delta;
}
assert !(Math.abs(error) >= Math.ulp(value)) : value; // Use ! for being tolerant to NaN.
return error;
}
/**
* Returns {@code true} if this {@code DoubleDouble} is equal to zero.
*
* @return {@code true} if this {@code DoubleDouble} is equal to zero.
*
* @see #ZERO
*/
public boolean isZero() {
return value == 0 && error == 0;
}
/**
* Returns {@code true} if this {@code DoubleDouble} is not a number.
*
* @return {@code true} if this {@code DoubleDouble} is not a number.
*
* @see #NaN
*/
public boolean isNaN() {
return Double.isNaN(value) || Double.isNaN(error);
}
/**
* Returns the sum of the given numbers,
* to be used only when {@code abs(a) >= abs(b)}.
*
* <p>Source: [Hida &amp; al.] page 4 algorithm 3, itself reproduced from [Shewchuk] page 312.</p>
*
* @param a the first number to add.
* @param b the second number to add, which must be smaller than {@code a}.
* @return sum of the given numbers.
*/
static DoubleDouble quickSum(final double a, final double b) {
final double value = a + b;
return new DoubleDouble(value, b - (value - a));
}
/**
* Returns the sum of the given numbers.
* The double-double accuracy is useful when the largest value is exact in base 2.
* A typical example is: 1 - (small value).
*
* <p>Source: [Hida &amp; al.] page 4 algorithm 4, itself reproduced from [Shewchuk] page 314.</p>
*
* @param a the first number to add.
* @param b the second number to add.
* @return sum of the given numbers.
*/
public static DoubleDouble sum(final double a, final double b) {
final double value = a + b;
final double v = value - a;
return new DoubleDouble(value, (a - (value - v)) + (b - v));
}
/**
* Returns the product of the given numbers.
* Note that unless the given arguments are exact in base 2,
* the result is not more accurate than a {@code double}.
*
* <p>Source: [Hida &amp; al.] page 5 algorithm 7, itself reproduced from [Shewchuk] page 326.
* This is the algorithm used when FMA instruction is available. For an algorithm without FMA,
* see [Hida &amp; al.] page 4 algorithm 6 (it is more complicated).</p>
*
* @param a the first number to multiply.
* @param b the second number to multiply.
* @return product of the given numbers.
*/
public static DoubleDouble product(final double a, final double b) {
final double value = a * b;
return new DoubleDouble(value, Math.fma(a, b, -value)); // Really needs the `fma(…)` precision here.
}
/**
* {@return {@code 1/this}}.
*/
public DoubleDouble inverse() {
return ONE.divide(this);
}
/**
* {@return {@code -this}}.
*/
public DoubleDouble negate() {
return new DoubleDouble(-value, -error);
}
/**
* Adds another double-double value to this {@code DoubleDouble}.
*
* @param other the other value to add to this {@code DoubleDouble}.
* @return the sum of {@code this} with the given number.
*/
public DoubleDouble add(final DoubleDouble other) {
return add(other.value, other.error);
}
/**
* Adds a {@code Number} value to this {@code DoubleDouble}. If the given number is an instance
* of {@code DoubleDouble} or {@link Fraction}, then the error term will be taken in account.
*
* @param other the other value to add to this {@code DoubleDouble}.
* @param decimal whether {@code float} and {@code double} values were intended to be exact in base 10.
* @return the sum of {@code this} with the given number.
*/
public DoubleDouble add(final Number other, final boolean decimal) {
return add(of(other, decimal));
}
/**
* Adds a {@code int} value to this {@code DoubleDouble}.
*
* @param other the other value to add to this {@code DoubleDouble}.
* @return the sum of {@code this} with the given number.
*/
public DoubleDouble add(final int other) {
return add(other, 0);
}
/**
* Adds a {@code long} value to this {@code DoubleDouble}.
*
* @param other the other value to add to this {@code DoubleDouble}.
* @return the sum of {@code this} with the given number.
*/
public DoubleDouble add(final long other) {
return add(of(other));
}
/**
* Adds a {@code double} value to this {@code DoubleDouble}.
* If {@code decimal} is {@code true}, then an error term is inferred for well-known values.
*
* @param other the other value to add to this {@code DoubleDouble}.
* @param decimal whether the value was intended to be exact in base 10.
* @return the sum of {@code this} with the given number.
*/
public DoubleDouble add(final double other, final boolean decimal) {
return add(other, decimal ? errorForWellKnownValue(other) : 0);
}
/**
* Adds another double-double value to this {@code DoubleDouble}.
*
* <h4>Implementation</h4>
* If <var>a</var> and <var>b</var> are {@code DoubleDouble} instances, then:
*
* <blockquote>(a + b)</blockquote>
*
* can be computed as:
*
* <blockquote>(a.value + a.error) + (b.value + b.error)<br>
* = (a.value + b.value) + (a.error + b.error)</blockquote>
*
* keeping in mind that the result of (a.value + b.value) has itself an error
* which needs to be added to (a.error + b.error). In Java code:
*
* {@snippet lang="java" :
* final double thisError = this.error;
* setToSum(value, otherValue);
* error += thisError;
* error += otherError;
* setToQuickSum(value, error);
* }
*
* @param otherValue the other value to add to this {@code DoubleDouble}.
* @param otherError the error of the other value to add to this {@code DoubleDouble}.
* @return the sum of {@code this} with the given number.
*/
private DoubleDouble add(final double otherValue, final double otherError) {
double s = value + otherValue;
double v = s - value;
double e = (value - (s - v)) + (otherValue - v) + (error + otherError);
if (s == 0 && e != 0) {
/*
* The two values almost cancelled, only their error terms are different.
* The number of significand bits (mantissa) in the IEEE `double` representation is 52,
* not counting the hidden bit. So estimate the accuracy of the double-double number as
* the accuracy of the `double` value (which is 1 ULP) scaled as if we had 52 additional
* significand bits (we ignore some more bits if ZERO_THRESHOLD is greater than 0).
* If the error is not greater than that value, then assume that it is not significant.
*/
if (Math.abs(e) <= Math.scalb(Math.ulp(otherValue), ZERO_THRESHOLD - (DOUBLE_PRECISION - 1))) {
return new DoubleDouble(s, 0);
}
}
return quickSum(s, e);
}
/**
* Subtracts another double-double value from this {@code DoubleDouble}.
*
* @param other the other value to subtract from this value.
* @return the difference between {@code this} and the given number.
*/
public DoubleDouble subtract(final DoubleDouble other) {
return add(-other.value, -other.error);
}
/**
* Subtracts a {@code Number} from this {@code DoubleDouble}. If the given number is an instance
* of {@code DoubleDouble} or {@link Fraction}, then the error term will be taken in account.
*
* @param other the other value to subtract from this {@code DoubleDouble}.
* @param decimal whether {@code float} and {@code double} values were intended to be exact in base 10.
* @return the difference between {@code this} and the given number.
*/
public DoubleDouble subtract(final Number other, final boolean decimal) {
return subtract(of(other, decimal));
}
/**
* Subtracts an {@code int} from this {@code DoubleDouble}.
*
* @param other the other value to subtract from this {@code DoubleDouble}.
* @return the difference between {@code this} and the given number.
*/
public DoubleDouble subtract(final int other) {
return add(-((double) other), 0);
}
/**
* Subtracts a {@code long} from this {@code DoubleDouble}.
*
* @param other the other value to subtract from this {@code DoubleDouble}.
* @return the difference between {@code this} and the given number.
*/
public DoubleDouble subtract(final long other) {
return subtract(of(other));
}
/**
* Subtracts a {@code double} from this {@code DoubleDouble} with a default error term.
* If {@code decimal} is {@code true}, then an error term is inferred for well-known values.
*
* @param other the other value to subtract from this {@code DoubleDouble}.
* @param decimal whether the value was intended to be exact in base 10.
* @return the difference between {@code this} and the given number.
*/
public DoubleDouble subtract(double other, final boolean decimal) {
other = -other;
return add(other, decimal ? errorForWellKnownValue(other) : 0);
}
/**
* Multiplies this {@code DoubleDouble} by another double-double value.
*
* @param other the other value to multiply by this value.
* @return the product of {@code this} with the given number.
*/
public DoubleDouble multiply(final DoubleDouble other) {
return multiply(other.value, other.error);
}
/**
* Multiplies this {@code DoubleDouble} by a {@code Number}. If the given number is an instance
* of {@code DoubleDouble} or {@link Fraction}, then the error term will be taken in account.
*
* @param other the other value to multiply by this {@code DoubleDouble}.
* @param decimal whether {@code float} and {@code double} values were intended to be exact in base 10.
* @return the product of {@code this} with the given number.
*/
public DoubleDouble multiply(final Number other, final boolean decimal) {
return multiply(of(other, decimal));
}
/**
* Multiplies this {@code DoubleDouble} by an {@code int}.
*
* @param other the other value to multiply by this {@code DoubleDouble}.
* @return the product of {@code this} with the given number.
*/
public DoubleDouble multiply(final int other) {
return multiply(other, 0);
}
/**
* Multiplies this {@code DoubleDouble} by a {@code long}.
*
* @param other the other value to multiply by this {@code DoubleDouble}.
* @return the product of {@code this} with the given number.
*/
public DoubleDouble multiply(final long other) {
return multiply(of(other));
}
/**
* Multiplies this {@code DoubleDouble} by a {@code double} with a default error term.
* If {@code decimal} is {@code true}, then an error term is inferred for well-known values.
*
* @param other the other value to multiply by this {@code DoubleDouble}.
* @param decimal whether the value was intended to be exact in base 10.
* @return the product of {@code this} with the given number.
*/
public DoubleDouble multiply(final double other, final boolean decimal) {
return multiply(other, decimal ? errorForWellKnownValue(other) : 0);
}
/**
* Multiplies this {@code DoubleDouble} by another double-double value.
*
* <h4>Implementation</h4>
* If <var>a</var> and <var>b</var> are {@code DoubleDouble} instances, then:
*
* <blockquote>(a * b)</blockquote>
*
* can be computed as:
*
* <blockquote>(a.value + a.error) * (b.value + b.error)<br>
* = (a.value * b.value) + (a.error * b.value) + (a.value * b.error) + (a.error * b.error)<br>
* ≅ (a.value * b.value) + (a.error * b.value) + (a.value * b.error)</blockquote>
*
* The first term is the main product. All other terms are added to the error, keeping in mind that the main
* product has itself an error. The last term (the product of errors) is ignored because presumed very small.
*
* @param otherValue the other value by which to multiply this {@code DoubleDouble}.
* @param otherError the error of the other value by which to multiply this {@code DoubleDouble}.
* @return the product of {@code this} with the given number.
*/
private DoubleDouble multiply(final double otherValue, final double otherError) {
double v = value * otherValue, e;
e = Math.fma(value, otherValue, -v); // Really needs the `fma(…)` precision here.
e = Math.fma(otherError, value, e);
e = Math.fma(otherValue, error, e);
return quickSum(v, e);
}
/**
* Divides this {@code DoubleDouble} by another double-double value.
*
* @param other the other value to by which to divide this value.
* @return the ratio between {@code this} and the given number.
*/
public DoubleDouble divide(final DoubleDouble other) {
return divide(other.value, other.error);
}
/**
* Divides this {@code DoubleDouble} by a {@code Number}. If the given number is an instance
* of {@code DoubleDouble} or {@link Fraction}, then the error term will be taken in account.
*
* @param other the other value by which to divide this {@code DoubleDouble}.
* @param decimal whether {@code float} and {@code double} values were intended to be exact in base 10.
* @return the ratio between {@code this} and the given number.
*/
public DoubleDouble divide(final Number other, final boolean decimal) {
return divide(of(other, decimal));
}
/**
* Divides this {@code DoubleDouble} by an {@code int}.
*
* @param other the other value by which to divide this {@code DoubleDouble}.
* @return the ratio between {@code this} and the given number.
*/
public DoubleDouble divide(final int other) {
return divide(other, 0);
}
/**
* Divides this {@code DoubleDouble} by a {@code long}.
*
* @param other the other value by which to divide this {@code DoubleDouble}.
* @return the ratio between {@code this} and the given number.
*/
public DoubleDouble divide(final long other) {
return divide(of(other));
}
/**
* Divides this {@code DoubleDouble} by a {@code double} with a default error term.
* If {@code decimal} is {@code true}, then an error term is inferred for well-known values.
*
* @param other the other value by which to divide this {@code DoubleDouble}.
* @param decimal whether the value was intended to be exact in base 10.
* @return the ratio between {@code this} and the given number.
*/
public DoubleDouble divide(final double other, final boolean decimal) {
return divide(other, decimal ? errorForWellKnownValue(other) : 0);
}
/**
* Divides this {@code DoubleDouble} by another double-double value.
*
* <h4>Implementation</h4>
* If <var>a</var> and <var>b</var> are {@code DoubleDouble} instances, then we estimate:
*
* <blockquote>(a / b) = (a.value / b.value) + remainder / b</blockquote>
*
* where:
*
* <blockquote>remainder = a - b * (a.value / b.value)</blockquote>
*
* @param otherValue the other value by which to divide this {@code DoubleDouble}.
* @param otherError the error of the other value by which to divide this {@code DoubleDouble}.
* @return the ratio between {@code this} and the given number.
*/
private DoubleDouble divide(final double otherValue, final double otherError) {
if (DISABLED) {
return new DoubleDouble(value / otherValue, 0);
}
/*
* The `b * (a.value / b.value)` part in the method javadoc.
* The result should be `a` ± some error to be determined.
*/
double pv, pe, s;
final double quotient = value / otherValue;
pe = Math.fma(quotient, otherValue, -value); // Really needs the `fma(…)` precision here.
pe = Math.fma(quotient, otherError, pe);
pv = value + pe; // `quickSum(value, r)` inlined on 3 lines.
s = value - pv;
pe += s;
/*
* Compute the remainder as `a - above_product` where the product is the (pv + pe) pair.
* Code below is a call to `sum(-pv, -pe)` inlined and without final `quickSum(…)` call.
*/
final double v = s - value;
double e = (value - (s - v)) - (pv + v) + (error - pe);
/*
* Adds the `remainder / b` term, using `remainder / b.value` as an approximation
* (otherwise we would have to invoke this method recursively). The approximation
* is assumed okay because the second term is small compared to the first one.
*/
return quickSum(quotient, (s + e) / otherValue);
}
/**
* Computes (1-x)/(1+x) where <var>x</var> is {@code this}.
* This pattern occurs in map projections.
*
* @return (1-x)/(1+x).
*/
public DoubleDouble ratio_1m_1p() {
return ONE.subtract(this).divide(add(1));
}
/**
* Returns {@code this} × 2ⁿ. Typical usages are
* {@code scalb(1)} for an efficient multiplication by 2 and
* {@code scalb(-1)} for an efficient division by 2.
*
* @param n power of 2 used to scale {@code this}.
* @return {@code this} × 2ⁿ.
*/
public DoubleDouble scalb(final int n) {
return new DoubleDouble(Math.scalb(value, n), Math.scalb(error, n));
}
/**
* Computes the square of this value.
*
* @return {@code this * this}.
*/
public DoubleDouble square() {
return multiply(value, error);
}
/**
* Sets this double-double value to its square root.
*
* <h4>Implementation</h4>
* This method searches for a {@code (r + ε)} value where:
*
* <blockquote>(r + ε)² = {@linkplain #value} + {@linkplain #error}</blockquote>
*
* If we could compute {@code r = sqrt(value + error)} with enough precision, then ε would be 0.
* But with the {@code double} type, we can only estimate {@code r ≈ sqrt(value)}. However, since
* that <var>r</var> value should be close to the "true" value, then ε should be small.
*
* <blockquote>value + error = (r + ε)² = r² + 2rε + ε²</blockquote>
*
* Neglecting ε² on the assumption that |ε| ≪ |r|:
*
* <blockquote>value + error ≈ r² + 2rε</blockquote>
*
* Isolating ε:
*
* <blockquote>ε ≈ (value + error - r²) / (2r)</blockquote>
*
* @return the square root of this value.
*/
public DoubleDouble sqrt() {
if (value == 0) {
return ZERO;
}
double r = Math.sqrt(value);
DoubleDouble t = product(r, r);
t = t.subtract(this);
t = t.divide(-2*r, 0); // Multiplication by 2 does not cause any precision lost.
return quickSum(r, t.value);
}
/**
* Computes c₀ + c₁x + c₂x² + c₃x³ + c₄x⁴ + … where <var>x</var> is {@code this}.
* The given <var>c</var> coefficients are presumed accurate in base 2
* (i.e. this method does not try to apply a correction for base 10).
*
* @param coefficients the {@code c} coefficients. The array length must be at least 1.
* @return the series sum.
*/
public DoubleDouble series(final double... coefficients) {
DoubleDouble sum = new DoubleDouble(coefficients[0], 0);
DoubleDouble xn = this;
for (int i=1; i < coefficients.length; i++) {
sum = sum.add(xn.multiply(coefficients[i], 0));
xn = multiply(xn);
}
return sum;
}
/**
* Compares this value with the given value for order.
*
* @param other the value to be compared.
*/
@Override
public int compareTo(final DoubleDouble other) {
int c = Double.compare(value, other.value);
if (c == 0) {
c = Double.compare(error, other.error);
}
return c;
}
/**
* Compares this number with the given object for equality.
*
* @param obj the other object to compare with this number.
* @return {@code true} if both object are equal.
*/
@Override
public boolean equals(Object obj) {
if (obj instanceof DoubleDouble) {
final DoubleDouble other = (DoubleDouble) obj;
return Numerics.equals(value, other.value) &&
Numerics.equals(error, other.error);
}
return false;
}
/**
* Returns a hash code value for this number.
*
* @return a hash code value.
*/
@Override
public int hashCode() {
return Long.hashCode(Double.doubleToLongBits(value) ^ Double.doubleToLongBits(error));
}
/**
* Returns a string representation of this number.
* This method prints only {@link #value} digits because the {@link #error} term is not always significant.
* For example SIS sometime uses {@code DoubleDouble} in calculations involving trigonometric operations
* where operands have only {@code double} precision. Because users may see {@code DoubleDouble} values
* returned by {@code MatrixSIS}, we want to avoid misleading them with non-realistic precision.
*
* @return a string representation of this number.
*/
@Override
public String toString() {
return Double.toString(doubleValue());
}
}