Accumulator.java

/*
 * Copyright (C) 2018 Alberto Irurueta Carro (alberto@irurueta.com)
 *
 * Licensed 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 com.irurueta.navigation.geodesic;

/**
 * An accumulator for sums.
 * This allows many double precision numbers to be added together with twice the normal
 * precision. Thus, the effective precision of the sum is 106 bits or about 32 decimal places.
 * The implementation follows J. R. Shewchuk,
 * <a href="https://doi.org/10.1007/PL00009321">Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric
 * Predicates</a>, Discrete &amp; Computational Geometry 18(3) 305&ndash;363 (1997).
 * In the documentation of the member functions, <i>sum</i> stands for the value currently held in the
 * accumulator.
 */
public class Accumulator {

    /**
     * s + t accumulators for the sum.
     */
    private double s;

    /**
     * s + t accumulators for the sum.
     */
    private double t;

    /**
     * Constructor from a double.
     *
     * @param y set <i>sum</i> = <i>y</i>.
     */
    public Accumulator(final double y) {
        s = y;
        t = 0;
    }

    /**
     * Constructor from another Accumulator.
     *
     * @param a set <i>sum</i> = <i>a</i>.
     */
    public Accumulator(final Accumulator a) {
        s = a.s;
        t = a.t;
    }

    /**
     * Sets the value to a double.
     *
     * @param y set <i>sum</i> = <i>y</i>.
     */
    public void set(final double y) {
        s = y;
        t = 0;
    }

    /**
     * Returns the value held in the accumulator.
     *
     * @return <i>sum</i>.
     */
    public double getSum() {
        return s;
    }

    /**
     * Returns the result of adding a number to <i>sum</i> (but don't change <i>sum</i>).
     *
     * @param y the number to be added to the sum.
     * @return <i>sum</i> + <i>y</i>.
     */
    public double sum(final double y) {
        final var a = new Accumulator(this);
        a.add(y);
        return a.s;
    }

    /**
     * Add a number to the accumulator.
     *
     * @param y set <i>sum</i> += <i>y</i>.
     */
    public void add(double y) {
        // Here's Shewchuk's solution...

        // hold exact sum as [s, t, u]
        double u;

        // accumulate starting at least significant end
        var r = GeoMath.sum(y, t);
        y = r.getFirst();
        u = r.getSecond();

        r = GeoMath.sum(y, s);
        s = r.getFirst();
        t = r.getSecond();

        // Start is mS, mT decreasing and non-adjacent. Sum is now (s + t + u) exactly with s, t, u
        // non-adjacent and in decreasing order (except for possible zeros). The following code tries
        // to normalize the result.
        // Ideally, we want mS = round(s + t + u) and mU = round(s + t + u - mS). The following does an
        // approximate job (and maintains the decreasing non-adjacent property). Here are two "failures"
        // using 3-bit floats:

        // Case 1: mS is not equal to round(s + t + u) -- off by 1 ulp
        // [12, -1] - 8 -> [4, 0, -1] -> [4, -1] = 3 should be [3, 0] = 3

        // Case 2: mS + mT is not as close to s + t + u as it should be
        // [64, 5] + 4 -> [64, 8, 1] -> [64, 8] = 72 (off by 1)
        //               should be [80, -7] = 73 (exact)

        // "Fixing" these problems is probably not worth the expense. The representation inevitably
        // leads to small errors in the accumulated values.
        // The additional errors illustrated here amount to 1 ulp of the less significant word during
        // each addition to the Accumulator and an additional possible error of 1 ulp in the reported sum.

        // Incidentally, the "ideal" representation described above is not canonical, because
        // mS = round(mS + mT) may not be true. For example with 3-bit floats:

        // [128, 16] + 1 -> [160, -16] -- 160 = round(145).
        // But [160, 0] - 16 -> [128, 16] -- 128 = round(144).

        if (s == 0) {
            // this implies t == 0, so result is u.
            s = u;
        } else {
            // otherwise just accumulate u to t.
            t += u;
        }
    }

    /**
     * Negate an accumulator.
     * Set <i>sum</i> = &minus;<i>sum</i>.
     */
    public void negate() {
        s = -s;
        t = -t;
    }
}