/*
 * 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.internal.util;

import java.util.Arrays;
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;


/**
 * Basic arithmetic methods for extended precision numbers using the <cite>double-double</cite> 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:
 *
 * {@preformat java
 *     BigDecimal decimal = new BigDecimal(dd.value).add(new BigDecimal(dd.error));
 * }
 *
 * <div class="section">Impact of availability of FMA instructions</div>
 * When allowed to use <cite>fused multiply-add</cite> (FMA) instruction added in JDK9
 * (see <a href="https://issues.apache.org/jira/browse/SIS-136">SIS-136</a> on Apache SIS JIRA),
 * then the following methods should be revisited:
 *
 * <ul>
 *   <li>{@link #setToProduct(double, double)} - revisit with [Hida &amp; al.] algorithm 7.</li>
 * </ul>
 *
 * @author  Martin Desruisseaux (Geomatys)
 * @version 1.0
 *
 * @see <a href="http://en.wikipedia.org/wiki/Double-double_%28arithmetic%29#Double-double_arithmetic">Wikipedia: Double-double arithmetic</a>
 *
 * @since 0.4
 * @module
 */
public final class DoubleDouble extends Number {
    /**
     * 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>Since 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. This is normal. 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 {@code double-double}
     * arithmetic if the results are to be compared for strict equality.</p>
     */
    public static final boolean DISABLED = false;

    /**
     * When computing <var>a</var> - <var>b</var> as a double-double (106 significand bits) value,
     * if the amount of non-zero significand bits is equals or lower than {@code ZERO_THRESHOLD+1},
     * consider the result as zero.
     */
    private static final int ZERO_THRESHOLD = 2;

    /**
     * The split constant used as part of multiplication algorithms. The split algorithm is as below
     * (we have to inline it in multiplication methods because Java can not return multi-values):
     *
     * {@preformat java
     *     private void split(double a) {
     *         double t   = SPLIT * a;
     *         double ahi = t - (t - a);
     *         double alo = a - ahi;
     *     }
     * }
     *
     * <p>Source: [Hida &amp; al.] page 4 algorithm 5, itself reproduced from [Shewchuk] page 325.</p>
     */
    private static final double SPLIT = (1 << 27) + 1;

    /**
     * Maximal value that can be handled by {@link #multiply(double, double)}.
     * If a multiplication is using a value greater than {@code MAX_VALUE},
     * then the result will be infinity or NaN.
     */
    public static final double MAX_VALUE = Double.MAX_VALUE / SPLIT;

    /**
     * Pre-defined constants frequently used in SIS, sorted in increasing order. This table contains only
     * constants that can not 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
    };

    /**
     * The main value, minus the {@link #error}.
     */
    public double value;

    /**
     * The error that shall be added to the main {@link #value} in order to get the
     * <cite>"real"</cite> (actually <cite>"the most accurate that we can"</cite>) value.
     */
    public double error;

    /**
     * Creates a new value initialized to zero.
     */
    public DoubleDouble() {
    }

    /**
     * Creates a new value initialized to the given value.
     *
     * @param  other  the other value to copy.
     */
    public DoubleDouble(final DoubleDouble other) {
        value = other.value;
        error = other.error;
    }

    /**
     * Creates a new value initialized to 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  otherValue  the initial value.
     */
    public DoubleDouble(Number otherValue) {
        if (otherValue instanceof Fraction) {
            value = ((Fraction) otherValue).denominator;
            inverseDivide(((Fraction) otherValue).numerator);
        } else {
            if (otherValue instanceof BigInteger) {
                otherValue = new BigDecimal((BigInteger) otherValue, MathContext.DECIMAL128);
            }
            value = otherValue.doubleValue();
            if (otherValue instanceof DoubleDouble) {
                error = ((DoubleDouble) otherValue).error;
            } else if (otherValue instanceof BigDecimal) {
                // Really need new BigDecimal(value) below, not BigDecimal.valueOf(value).
                error = ((BigDecimal) otherValue).subtract(new BigDecimal(value), MathContext.DECIMAL64).doubleValue();
            } else {
                error = errorForWellKnownValue(value);
            }
        }
    }

    /**
     * Returns {@code true} if the given value is one of the special cases recognized by the
     * {@link DoubleDouble(Number)} constructor. Those special cases should rarely occur, so
     * we do not complicate the code with optimized code paths.
     *
     * @param  value  the value to test.
     * @return {@code true} if it is worth to convert the given value to a {@code DoubleDouble}.
     *
     * @since 0.8
     */
    public static boolean shouldConvert(final Number value) {
        return (value instanceof Fraction) || (value instanceof BigInteger) || (value instanceof BigDecimal);
    }

    /**
     * Creates a new instance initialized to the given long integer.
     *
     * @param  value  the long integer value to wrap.
     */
    public DoubleDouble(final long value) {
        this.value = value;
        this.error = (value - (long) this.value);
    }

    /**
     * Creates a new instance initialized to the given value verbatim, without inferring an error term for double-double arithmetic.
     * We use this constructor 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. It should also be used when the value is known to have
     * an exact representation as a {@code double} primitive type.
     *
     * @param  value  the value to wrap in a {@code DoubleDouble} instance.
     *
     * @see #createAndGuessError(double)
     */
    public DoubleDouble(final double value) {
        this.value = value;
    }

    /**
     * Creates a new value initialized to the given value and error.
     * It is caller's responsibility to ensure that the (value, error) pair is normalized.
     *
     * @param  value  the initial value.
     * @param  error  the initial error.
     */
    public DoubleDouble(final double value, final double error) {
        this.value = value;
        this.error = error;
        assert !(Math.abs(error) >= Math.ulp(value)) : this;            // Use ! for being tolerant to NaN.
    }

    /**
     * Returns the given value as a {@code DoubleDouble}. This method returns the given instance
     * directly if it can be safely casted to {@code DoubleDouble}.
     *
     * @param  value  the value to cast or to copy, or {@code null}.
     * @return the value as a {@code DoubleDouble} (may be the same instance than the given argument),
     *         or {@code null} if the given value was null.
     *
     * @since 0.8
     */
    public static DoubleDouble castOrCopy(final Number value) {
        return (value == null || value instanceof DoubleDouble) ? (DoubleDouble) value : new DoubleDouble(value);
    }

    /**
     * Creates a new value initialized to the given value and an error term inferred by
     * {@link #errorForWellKnownValue(double)}.
     *
     * <b>Tip:</b> if the other value is known to be an integer or a power of 2, then invoking
     * <code>{@linkplain #DoubleDouble(double) DoubleDouble}(value)</code> is more efficient.
     *
     * @param  value  the initial value.
     * @return an instance initialized to the given value and a default error term.
     */
    public static DoubleDouble createAndGuessError(final double value) {
        return new DoubleDouble(value, errorForWellKnownValue(value));
    }

    /**
     * Returns a new {@code DoubleDouble} instance initialized to the conversion factor
     * from radians to angular degrees.
     *
     * @return an instance initialized to the 57.2957795130823208767981548141052 value.
     */
    public static DoubleDouble createRadiansToDegrees() {
        return new DoubleDouble(57.2957795130823208767981548141052, -1.9878495670576283E-15);
    }

    /**
     * Returns a new {@code DoubleDouble} instance initialized to the conversion factor
     * from angular degrees to radians.
     *
     * @return an instance initialized to the 0.01745329251994329576923690768488613 value.
     */
    public static DoubleDouble createDegreesToRadians() {
        return new DoubleDouble(0.01745329251994329576923690768488613, 2.9486522708701687E-19);
    }

    /**
     * Returns a new {@code DoubleDouble} instance initialized to the conversion factor
     * from arc-seconds to radians.
     *
     * @return an instance initialized to the 0.000004848136811095359935899141023579480 value.
     */
    public static DoubleDouble createSecondsToRadians() {
        return new DoubleDouble(0.000004848136811095359935899141023579480, 9.320078015422868E-23);
    }

    /** @return {@link #value} + {@link #error}. */
    @Override public double doubleValue() {return value + error;}
    @Override public float  floatValue()  {return (float) doubleValue();}
    @Override public long   longValue()   {return Math.round(value) + (long) error;}
    @Override public int    intValue()    {return Math.toIntExact(longValue());}

    /**
     * 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)}.
     *
     * <div class="section">Rational</div>
     * 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 can not 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 later case,
     *         the base 2 representation of the given value is assumed to be accurate enough.
     */
    public static double errorForWellKnownValue(final double value) {
        if (DISABLED) return 0;
        final int i = Arrays.binarySearch(VALUES, Math.abs(value));
        if (i >= 0) {
            return MathFunctions.xorSign(ERRORS[i], value);
        }
        final double delta = DecimalFunctions.deltaForDoubleToDecimal(value);
        return Double.isNaN(delta) ? 0 : delta;
    }

    /**
     * Returns {@code true} if this {@code DoubleDouble} is equals to zero.
     *
     * @return {@code true} if this {@code DoubleDouble} is equals to zero.
     */
    public boolean isZero() {
        return value == 0 && error == 0;
    }

    /**
     * Resets the {@link #value} and {@link #error} terms to zero.
     */
    public void clear() {
        value = 0;
        error = 0;
    }

    /**
     * Sets this {@code DoubleDouble} to the given 64-bits signed integer.
     *
     * @param  n  the value to set.
     */
    public void set(final long n) {
        value = n;
        error = n - (long) value;
    }

    /**
     * Sets this {@code DoubleDouble} to the same value than the given instance.
     *
     * @param  other  the instance to copy.
     */
    public void setFrom(final DoubleDouble other) {
        value = other.value;
        error = other.error;
    }

    /**
     * Sets the {@link #value} and {@link #error} terms to values read from the given array.
     * This is a convenience method for a frequently used operation, implemented as below:
     *
     * {@preformat java
     *   value = array[index];
     *   error = array[index + errorOffset];
     * }
     *
     * @param  array        the array from which to get the value and error.
     * @param  index        index of the value in the given array.
     * @param  errorOffset  offset to add to {@code index} in order to get the index of the error in the given array.
     */
    public void setFrom(final double[] array, final int index, final int errorOffset) {
        value = array[index];
        error = array[index + errorOffset];
    }

    /**
     * Equivalent to a call to {@code setToQuickSum(value, error)} inlined.
     * This is invoked after addition or multiplication operations.
     */
    final void normalize() {
        error += (value - (value += error));
        if (DISABLED) error = 0;
    }

    /**
     * Sets this {@code DoubleDouble} to 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}.
     */
    final void setToQuickSum(final double a, final double b) {
        value = a + b;
        error = b - (value - a);
        if (DISABLED) error = 0;
    }

    /**
     * Sets this {@code DoubleDouble} to the sum of the given numbers.
     *
     * <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.
     */
    public void setToSum(final double a, final double b) {
        value = a + b;
        final double v = value - a;
        error = (a - (value - v)) + (b - v);
        if (DISABLED) error = 0;
    }

    /**
     * Sets this {@code DoubleDouble} to the product of the given numbers.
     * The given numbers shall not be greater than {@value #MAX_VALUE} in magnitude.
     *
     * <p>Source: [Hida &amp; al.] page 4 algorithm 6, itself reproduced from [Shewchuk] page 326.</p>
     *
     * @param  a  the first number to multiply.
     * @param  b  the second number to multiply.
     */
    public void setToProduct(final double a, final double b) {
        value = a * b;
        double t = SPLIT * a;
        final double ahi = t - (t - a);
        final double alo = a - ahi;
        t = SPLIT * b;
        final double bhi = t - (t - b);
        final double blo = b - bhi;
        error = ((ahi*bhi - value) + ahi*blo + alo*bhi) + alo*blo;
        if (DISABLED) error = 0;
    }

    /**
     * Stores the {@link #value} and {@link #error} terms in the given array.
     * This is a convenience method for a frequently used operation, implemented as below:
     *
     * {@preformat java
     *   array[index] = value;
     *   array[index + errorOffset] = error;
     * }
     *
     * @param  array        the array where to store the value and error.
     * @param  index        index of the value in the given array.
     * @param  errorOffset  offset to add to {@code index} in order to get the index of the error in the given array.
     */
    public void storeTo(final double[] array, final int index, final int errorOffset) {
        array[index] = value;
        array[index + errorOffset] = error;
    }

    /**
     * Swaps two double-double values in the given array.
     *
     * @param  array        the array where to swap the values and errors.
     * @param  i0           index of the first value to swap.
     * @param  i1           index of the second value to swap.
     * @param  errorOffset  offset to add to the indices in order to get the error indices in the given array.
     *
     * @see org.apache.sis.util.ArraysExt#swap(double[], int, int)
     */
    public static void swap(final double[] array, int i0, int i1, final int errorOffset) {
        double t = array[i0];
        array[i0] = array[i1];
        array[i1] = t;
        t = array[i0 += errorOffset];
        array[i0] = array[i1 += errorOffset];
        array[i1] = t;
    }

    /**
     * Sets this number to {@code -this}.
     */
    public void negate() {
        value = -value;
        error = -error;
    }

    /**
     * Adds an other double-double value to this {@code DoubleDouble}.
     * This is a convenience method for:
     *
     * {@preformat java
     *    add(other.value, other.error);
     * }
     *
     * @param  other  the other value to add to this {@code DoubleDouble}.
     */
    public void add(final DoubleDouble other) {
        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}.
     */
    public void addGuessError(final Number other) {
        if (other instanceof DoubleDouble) {
            add((DoubleDouble) other);
        } else if (shouldConvert(other)) {
            add(new DoubleDouble(other));
        } else {
            addGuessError(other.doubleValue());
        }
    }

    /**
     * Adds a {@code double} value to this {@code DoubleDouble} with a default error term.
     * This is a convenience method for:
     *
     * {@preformat java
     *    add(otherValue, errorForWellKnownValue(otherValue));
     * }
     *
     * <b>Tip:</b> if the other value is known to be an integer or a power of 2, then invoking
     * <code>{@linkplain #add(double) add}(otherValue)</code> is more efficient.
     *
     * @param  otherValue  the other value to add to this {@code DoubleDouble}.
     */
    public void addGuessError(final double otherValue) {
        add(otherValue, errorForWellKnownValue(otherValue));
    }

    /**
     * Adds the given {@code double} value using
     * <a href="https://en.wikipedia.org/wiki/Kahan_summation_algorithm">Kahan summation algorithm</a>.
     * This can be used when {@code otherValue} is known to be smaller than {@link #value}.
     *
     * @param  y  the other value to add to this {@code DoubleDouble}.
     */
    public void addKahan(double y) {
        y += error;
        error = y + (value - (value += y));
    }

    /**
     * Adds a {@code double} value to this {@code DoubleDouble} without error term.
     * This is a convenience method for:
     *
     * {@preformat java
     *    add(otherValue, 0);
     * }
     *
     * @param  otherValue  the other value to add to this {@code DoubleDouble}.
     */
    public void add(final double otherValue) {
        add(otherValue, 0);
    }

    /**
     * Adds an other double-double value to this {@code DoubleDouble}.
     * The result is stored in this instance.
     *
     * <div class="section">Implementation</div>
     * 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:
     *
     * {@preformat 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}.
     */
    public void add(final double otherValue, final double otherError) {
        // Inline expansion of the code in above javadoc.
        double v = value;
        value += otherValue;
        error += v - (value + (v -= value)) + (otherValue + v);
        error += otherError;
        if (value == 0 && error != 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(error) <= Math.scalb(Math.ulp(otherValue), ZERO_THRESHOLD - Numerics.SIGNIFICAND_SIZE)) {
                error = 0;
                return;
            }
        }
        normalize();
    }

    /**
     * Adds an other double-double value to this {@code DoubleDouble}, reading the values from an array.
     * This is a convenience method for a frequently used operation, implemented as below:
     *
     * {@preformat java
     *    add(array[index], array[index + errorOffset]);
     * }
     *
     * @param  array        the array from which to get the value and error.
     * @param  index        index of the value in the given array.
     * @param  errorOffset  offset to add to {@code index} in order to get the index of the error in the given array.
     */
    public void add(final double[] array, final int index, final int errorOffset) {
        add(array[index], array[index + errorOffset]);
    }

    /**
     * Subtracts an other double-double value from this {@code DoubleDouble}.
     * This is a convenience method for:
     *
     * {@preformat java
     *    subtract(other.value, other.error);
     * }
     *
     * @param  other  the other value to subtract from this value.
     */
    public void subtract(final DoubleDouble other) {
        subtract(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}.
     */
    public void subtractGuessError(final Number other) {
        if (other instanceof DoubleDouble) {
            subtract((DoubleDouble) other);
        } else if (shouldConvert(other)) {
            subtract(new DoubleDouble(other));
        } else {
            subtractGuessError(other.doubleValue());
        }
    }

    /**
     * Subtracts a {@code double} from this {@code DoubleDouble} with a default error term.
     * This is a convenience method for:
     *
     * {@preformat java
     *    subtract(otherValue, errorForWellKnownValue(otherValue));
     * }
     *
     * <b>Tip:</b> if the other value is known to be an integer or a power of 2, then invoking
     * <code>{@linkplain #subtract(double) subtract}(otherValue)</code> is more efficient.
     *
     * @param  otherValue  the other value to subtract from this {@code DoubleDouble}.
     */
    public void subtractGuessError(final double otherValue) {
        subtract(otherValue, errorForWellKnownValue(otherValue));
    }

    /**
     * Subtracts a {@code double} from this {@code DoubleDouble} without error term.
     * This is a convenience method for:
     *
     * {@preformat java
     *    subtract(otherValue, 0);
     * }
     *
     * @param  otherValue  the other value to subtract from this {@code DoubleDouble}.
     */
    public void subtract(final double otherValue) {
        subtract(otherValue, 0);
    }

    /**
     * Subtracts an other double-double value from this {@code DoubleDouble}.
     * The result is stored in this instance.
     *
     * @param  otherValue  the other value to subtract from this {@code DoubleDouble}.
     * @param  otherError  the error of the other value to subtract from this {@code DoubleDouble}.
     */
    public void subtract(final double otherValue, final double otherError) {
        add(-otherValue, -otherError);
    }

    /**
     * Subtracts an other double-double value from this {@code DoubleDouble}, reading the values from an array.
     * This is a convenience method for a frequently used operation, implemented as below:
     *
     * {@preformat java
     *    subtract(array[index], array[index + errorOffset]);
     * }
     *
     * @param  array        the array from which to get the value and error.
     * @param  index        index of the value in the given array.
     * @param  errorOffset  offset to add to {@code index} in order to get the index of the error in the given array.
     */
    public void subtract(final double[] array, final int index, final int errorOffset) {
        subtract(array[index], array[index + errorOffset]);
    }

    /**
     * Multiplies this {@code DoubleDouble} by an other double-double value.
     * This is a convenience method for:
     *
     * {@preformat java
     *    multiply(other.value, other.error);
     * }
     *
     * @param  other  the other value to multiply by this value.
     */
    public void multiply(final DoubleDouble other) {
        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}.
     */
    public void multiplyGuessError(final Number other) {
        if (other instanceof DoubleDouble) {
            multiply((DoubleDouble) other);
        } else if (shouldConvert(other)) {
            multiply(new DoubleDouble(other));
        } else {
            multiplyGuessError(other.doubleValue());
        }
    }

    /**
     * Multiplies this {@code DoubleDouble} by a {@code double} with a default error term.
     * This is a convenience method for:
     *
     * {@preformat java
     *    multiply(otherValue, errorForWellKnownValue(otherValue));
     * }
     *
     * <b>Tip:</b> if the other value is known to be an integer or a power of 2, then invoking
     * <code>{@linkplain #multiply(double) multiply}(otherValue)</code> is more efficient.
     *
     * @param  otherValue  the other value to multiply by this {@code DoubleDouble}.
     */
    public void multiplyGuessError(final double otherValue) {
        multiply(otherValue, errorForWellKnownValue(otherValue));
    }

    /**
     * Multiplies this {@code DoubleDouble} by a {@code double} without error term.
     * This is a convenience method for:
     *
     * {@preformat java
     *    multiply(otherValue, 0);
     * }
     *
     * @param  otherValue  the other value to multiply by this {@code DoubleDouble}.
     */
    public void multiply(final double otherValue) {
        multiply(otherValue, 0);
    }

    /**
     * Multiplies this {@code DoubleDouble} by an other double-double value.
     * The result is stored in this instance.
     *
     * <div class="section">Implementation</div>
     * 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.
     * In Java code:
     *
     * {@preformat java
     *   final double thisValue = this.value;
     *   final double thisError = this.error;
     *   setToProduct(thisValue, otherValue);
     *   error += otherError * thisValue;
     *   error += otherValue * thisError;
     *   setToQuickSum(value, error);
     * }
     *
     * @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}.
     */
    public void multiply(final double otherValue, final double otherError) {
        final double thisValue = this.value;
        final double thisError = this.error;
        setToProduct(thisValue, otherValue);
        error += otherError * thisValue;
        error += otherValue * thisError;
        normalize();
    }

    /**
     * Multiplies this {@code DoubleDouble} by an other double-double value stored in the given array.
     * This is a convenience method for a frequently used operation, implemented as below:
     *
     * {@preformat java
     *    multiply(array[index], array[index + errorOffset]);
     * }
     *
     * @param  array        the array from which to get the value and error.
     * @param  index        index of the value in the given array.
     * @param  errorOffset  offset to add to {@code index} in order to get the index of the error in the given array.
     */
    public void multiply(final double[] array, final int index, final int errorOffset) {
        multiply(array[index], array[index + errorOffset]);
    }

    /**
     * Divides this {@code DoubleDouble} by an other double-double value.
     * This is a convenience method for:
     *
     * {@preformat java
     *    divide(other.value, other.error);
     * }
     *
     * @param  other  the other value to by which to divide this value.
     */
    public void divide(final DoubleDouble other) {
        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}.
     */
    public void divideGuessError(final Number other) {
        if (other instanceof DoubleDouble) {
            divide((DoubleDouble) other);
        } else if (shouldConvert(other)) {
            divide(new DoubleDouble(other));
        } else {
            divideGuessError(other.doubleValue());
        }
    }

    /**
     * Divides this {@code DoubleDouble} by a {@code double} with a default error term.
     * This is a convenience method for:
     *
     * {@preformat java
     *    divide(otherValue, errorForWellKnownValue(otherValue));
     * }
     *
     * <b>Tip:</b> if the other value is known to be an integer or a power of 2, then invoking
     * <code>{@linkplain #divide(double) divide}(otherValue)</code> is more efficient.
     *
     * @param  otherValue  the other value by which to divide this {@code DoubleDouble}.
     */
    public void divideGuessError(final double otherValue) {
        divide(otherValue, errorForWellKnownValue(otherValue));
    }

    /**
     * Divides this {@code DoubleDouble} by a {@code double} without error term.
     * This is a convenience method for:
     *
     * {@preformat java
     *    divide(otherValue, 0);
     * }
     *
     * @param  otherValue  the other value by which to divide this {@code DoubleDouble}.
     */
    public void divide(final double otherValue) {
        divide(otherValue, 0);
    }

    /**
     * Divides this {@code DoubleDouble} by an other double-double value.
     * The result is stored in this instance.
     *
     * @param  denominatorValue  the other value by which to divide this {@code DoubleDouble}.
     * @param  denominatorError  the error of the other value by which to divide this {@code DoubleDouble}.
     */
    public void divide(final double denominatorValue, final double denominatorError) {
        if (DISABLED) {
            value /= denominatorValue;
            error  = 0;
            return;
        }
        final double numeratorValue = value;
        final double numeratorError = error;
        value = denominatorValue;
        error = denominatorError;
        inverseDivide(numeratorValue, numeratorError);
    }

    /**
     * Divides this {@code DoubleDouble} by an other double-double value stored in the given array.
     * This is a convenience method for a frequently used operation, implemented as below:
     *
     * {@preformat java
     *    divide(array[index], array[index + errorOffset]);
     * }
     *
     * @param  array        the array from which to get the value and error.
     * @param  index        index of the value in the given array.
     * @param  errorOffset  offset to add to {@code index} in order to get the index of the error in the given array.
     */
    public void divide(final double[] array, final int index, final int errorOffset) {
        divide(array[index], array[index + errorOffset]);
    }

    /**
     * Divides the given double-double value by this {@code DoubleDouble}.
     * This is a convenience method for:
     *
     * {@preformat java
     *    inverseDivide(other.value, other.error);
     * }
     *
     * @param  other  the other value to divide by this value.
     */
    public void inverseDivide(final DoubleDouble other) {
        inverseDivide(other.value, other.error);
    }

    /**
     * Divides the given {@code Number} value by 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 divide by this {@code DoubleDouble}.
     */
    public void inverseDivideGuessError(final Number other) {
        if (other instanceof DoubleDouble) {
            inverseDivide((DoubleDouble) other);
        } else if (shouldConvert(other)) {
            inverseDivide(new DoubleDouble(other));
        } else {
            inverseDivideGuessError(other.doubleValue());
        }
    }

    /**
     * Divides the given {@code double} value by this {@code DoubleDouble} with a default error term.
     * This is a convenience method for:
     *
     * {@preformat java
     *    inverseDivide(numeratorValue, errorForWellKnownValue(numeratorValue));
     * }
     *
     * <b>Tip:</b> if the other value is known to be an integer or a power of 2, then invoking
     * <code>{@linkplain #inverseDivide(double) inverseDivide}(otherValue)</code> is more efficient.
     *
     * @param  numeratorValue  the other value to divide by this {@code DoubleDouble}.
     */
    public void inverseDivideGuessError(final double numeratorValue) {
        inverseDivide(numeratorValue, errorForWellKnownValue(numeratorValue));
    }

    /**
     * Divides the given {@code double} value by this {@code DoubleDouble} without error term.
     * This is a convenience method for:
     *
     * {@preformat java
     *    inverseDivide(numeratorValue, 0);
     * }
     *
     * @param  numeratorValue  the other value to divide by this {@code DoubleDouble}.
     */
    public void inverseDivide(final double numeratorValue) {
        inverseDivide(numeratorValue, 0);
    }

    /**
     * Divides the given double-double value by this {@code DoubleDouble}.
     * The result is stored in this instance.
     *
     * <div class="section">Implementation</div>
     * 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  numeratorValue  the other value to divide by this {@code DoubleDouble}.
     * @param  numeratorError  the error of the other value to divide by this {@code DoubleDouble}.
     */
    public void inverseDivide(final double numeratorValue, final double numeratorError) {
        if (DISABLED) {
            value = numeratorValue / value;
            error = 0;
            return;
        }
        final double denominatorValue = value;
        /*
         * The 'b * (a.value / b.value)' part in the method javadoc.
         */
        final double quotient = numeratorValue / denominatorValue;
        multiply(quotient);
        /*
         * Compute 'remainder' as 'a - above_product'.
         */
        final double productError = error;
        setToSum(numeratorValue, -value);
        error -= productError;  // Complete the above subtraction
        error += numeratorError;
        /*
         * 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 since the second term is small compared to the first one.
         */
        setToQuickSum(quotient, (value + error) / denominatorValue);
    }

    /**
     * Divides the given double-double value by this {@code DoubleDouble}.
     * This is a convenience method for a frequently used operation, implemented as below:
     *
     * {@preformat java
     *    inverseDivide(array[index], array[index + errorOffset]);
     * }
     *
     * @param  array        the array from which to get the value and error.
     * @param  index        index of the value in the given array.
     * @param  errorOffset  offset to add to {@code index} in order to get the index of the error in the given array.
     */
    public void inverseDivide(final double[] array, final int index, final int errorOffset) {
        inverseDivide(array[index], array[index + errorOffset]);
    }

    /**
     * Computes (1-x)/(1+x) where <var>x</var> is {@code this}.
     * This pattern occurs in map projections.
     */
    public void ratio_1m_1p() {
        final DoubleDouble numerator = new DoubleDouble(1d);
        numerator.subtract(this);
        add(1);
        inverseDivide(numerator);
    }

    /**
     * Computes the square of this value.
     */
    public void square() {
        multiply(value, error);
    }

    /**
     * Sets this double-double value to its square root.
     *
     * <div class="section">Implementation</div>
     * 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>
     */
    public void sqrt() {
        if (value != 0) {
            final double thisValue = this.value;
            final double thisError = this.error;
            double r = Math.sqrt(thisValue);
            setToProduct(r, r);
            subtract(thisValue, thisError);
            divide(-2*r);                           // Multiplication by 2 does not cause any precision lost.
            setToQuickSum(r, 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.
     */
    public void series(final double... coefficients) {
        final DoubleDouble x = new DoubleDouble(this);
        value = coefficients[0];
        error = 0;
        final int last = coefficients.length - 1;
        if (last >= 1) {
            final DoubleDouble xn = new DoubleDouble(x);
            final DoubleDouble t = new DoubleDouble(xn);
            for (int i=1; i<last; i++) {
                t.multiply(coefficients[i]);
                add(t);
                xn.multiply(x);
                t.setFrom(xn);
            }
            t.multiply(coefficients[last]);
            add(t);
        }
    }

    /**
     * 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));
    }

    /**
     * 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 string representation of this number for debugging purpose.
     * The returned string does not need to contains all digits that this {@code DoubleDouble} can handle.
     *
     * @return a string representation of this number.
     */
    @Override
    public String toString() {
        return new BigDecimal(value).add(new BigDecimal(error), MathContext.DECIMAL128).toString();
    }
}