Norm.java

/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package org.apache.commons.numbers.core;

/**
 * <a href="https://en.wikipedia.org/wiki/Norm_(mathematics)">Norm</a> functions.
 *
 * <p>The implementations provide increased numerical accuracy.
 * Algorithms primary source is the 2005 paper
 * <a href="https://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.2.1547">
 * Accurate Sum and Dot Product</a> by Takeshi Ogita, Siegfried M. Rump,
 * and Shin'ichi Oishi published in <em>SIAM J. Sci. Comput</em>.
 */
public enum Norm {
    /**
     * <a href="https://en.wikipedia.org/wiki/Norm_(mathematics)#Taxicab_norm_or_Manhattan_norm">
     *  Manhattan norm</a> (sum of the absolute values of the arguments).
     */
    L1(Norm::manhattan, Norm::manhattan, Norm::manhattan),
    /** Alias for {@link #L1}. */
    MANHATTAN(L1),
    /** <a href="https://en.wikipedia.org/wiki/Norm_(mathematics)#Euclidean_norm">Euclidean norm</a>. */
    L2(Norm::euclidean, Norm::euclidean, Norm::euclidean),
    /** Alias for {@link #L2}. */
    EUCLIDEAN(L2),
    /**
     * <a href="https://en.wikipedia.org/wiki/Norm_(mathematics)#Maximum_norm_(special_case_of:_infinity_norm,_uniform_norm,_or_supremum_norm)">
     *  Maximum norm</a> (maximum of the absolute values of the arguments).
     */
    LINF(Norm::maximum, Norm::maximum, Norm::maximum),
    /** Alias for {@link #LINF}. */
    MAXIMUM(LINF);

    /**
     * Threshold for scaling small numbers. This value is chosen such that doubles
     * set to this value can be squared without underflow. Values less than this must
     * be scaled up.
     */
    private static final double SMALL_THRESH = 0x1.0p-511;
    /**
     * Threshold for scaling large numbers. This value is chosen such that 2^31 doubles
     * set to this value can be squared and added without overflow. Values greater than
     * this must be scaled down.
     */
    private static final double LARGE_THRESH = 0x1.0p+496;
    /**
     * Threshold for scaling up a single value by {@link #SCALE_UP} without risking
     * overflow when the value is squared.
     */
    private static final double SAFE_SCALE_UP_THRESH = 0x1.0p-100;
    /** Value used to scale down large numbers. */
    private static final double SCALE_DOWN = 0x1.0p-600;
    /** Value used to scale up small numbers. */
    private static final double SCALE_UP = 0x1.0p+600;

    /** Threshold for the difference between the exponents of two Euclidean 2D input values
     * where the larger value dominates the calculation.
     */
    private static final int EXP_DIFF_THRESHOLD_2D = 54;

    /** Function of 2 arguments. */
    private final Two two;
    /** Function of 3 arguments. */
    private final Three three;
    /** Function of array argument. */
    private final Array array;

    /** Function of 2 arguments. */
    @FunctionalInterface
    private interface Two {
        /**
         * @param x Argument.
         * @param y Argument.
         * @return the norm.
         */
        double of(double x, double y);
    }
    /** Function of 3 arguments. */
    @FunctionalInterface
    private interface Three {
        /**
         * @param x Argument.
         * @param y Argument.
         * @param z Argument.
         * @return the norm.
         */
        double of(double x, double y, double z);
    }
    /** Function of array argument. */
    @FunctionalInterface
    private interface Array {
        /**
         * @param v Array of arguments.
         * @return the norm.
         */
        double of(double[] v);
    }

    /**
     * @param two Function of 2 arguments.
     * @param three Function of 3 arguments.
     * @param array Function of array argument.
     */
    Norm(Two two,
         Three three,
         Array array) {
        this.two = two;
        this.three = three;
        this.array = array;
    }

    /**
     * @param alias Alternative name.
     */
    Norm(Norm alias) {
        this.two = alias.two;
        this.three = alias.three;
        this.array = alias.array;
    }

    /**
     * Computes the norm.
     *
     * <p>Special cases:
     * <ul>
     *  <li>If either value is {@link Double#NaN}, then the result is {@link Double#NaN}.</li>
     *  <li>If either value is infinite and the other value is not {@link Double#NaN}, then
     *   the result is {@link Double#POSITIVE_INFINITY}.</li>
     * </ul>
     *
     * @param x Argument.
     * @param y Argument.
     * @return the norm.
     */
    public final double of(double x,
                           double y) {
        return two.of(x, y);
    }

    /**
     * Computes the norm.
     *
     * <p>Special cases:
     * <ul>
     *  <li>If any value is {@link Double#NaN}, then the result is {@link Double#NaN}.</li>
     *  <li>If any value is infinite and no value is not {@link Double#NaN}, then the
     *   result is {@link Double#POSITIVE_INFINITY}.</li>
     * </ul>
     *
     * @param x Argument.
     * @param y Argument.
     * @param z Argument.
     * @return the norm.
     */
    public final double of(double x,
                           double y,
                           double z) {
        return three.of(x, y, z);
    }

    /**
     * Computes the norm.
     *
     * <p>Special cases:
     * <ul>
     *  <li>If any value is {@link Double#NaN}, then the result is {@link Double#NaN}.</li>
     *  <li>If any value is infinite and no value is not {@link Double#NaN}, then the
     *   result is {@link Double#POSITIVE_INFINITY}.</li>
     * </ul>
     *
     * @param v Argument.
     * @return the norm.
     * @throws IllegalArgumentException if the array is empty.
     */
    public final double of(double[] v) {
        ensureNonEmpty(v);
        return array.of(v);
    }

    /** Computes the Manhattan norm.
     *
     * @param x first input value
     * @param y second input value
     * @return \(|x| + |y|\).
     *
     * @see #L1
     * @see #MANHATTAN
     * @see #of(double,double)
     */
    private static double manhattan(final double x,
                                    final double y) {
        return Math.abs(x) + Math.abs(y);
    }

    /** Computes the Manhattan norm.
     *
     * @param x first input value
     * @param y second input value
     * @param z third input value
     * @return \(|x| + |y| + |z|\)
     *
     * @see #L1
     * @see #MANHATTAN
     * @see #of(double,double,double)
     */
    private static double manhattan(final double x,
                                    final double y,
                                    final double z) {
        return Sum.of(Math.abs(x))
            .add(Math.abs(y))
            .add(Math.abs(z))
            .getAsDouble();
    }

    /** Computes the Manhattan norm.
     *
     * @param v input values
     * @return \(|v_0| + ... + |v_i|\)
     *
     * @see #L1
     * @see #MANHATTAN
     * @see #of(double[])
     */
    private static double manhattan(final double[] v) {
        final Sum sum = Sum.create();

        for (final double d : v) {
            sum.add(Math.abs(d));
        }

        return sum.getAsDouble();
    }

    /** Computes the Euclidean norm.
     * This implementation handles possible overflow or underflow.
     *
     * <p><strong>Comparison with Math.hypot()</strong>
     * While not guaranteed to return the same result, this method provides
     * similar error bounds as {@link Math#hypot(double, double)} (and may run faster on
     * some JVM).
     *
     * @param x first input
     * @param y second input
     * @return \(\sqrt{x^2 + y^2}\).
     *
     * @see #L2
     * @see #EUCLIDEAN
     * @see #of(double,double)
     */
    private static double euclidean(final double x,
                                    final double y) {
        final double xabs = Math.abs(x);
        final double yabs = Math.abs(y);

        final double max;
        final double min;
        // the compare method considers NaN greater than other values, meaning that our
        // check for if the max is finite later on will detect NaNs correctly
        if (Double.compare(xabs, yabs) > 0) {
            max = xabs;
            min = yabs;
        } else {
            max = yabs;
            min = xabs;
        }

        // if the max is not finite, then one of the inputs must not have
        // been finite
        if (!Double.isFinite(max)) {
            // let the standard multiply operation determine whether to return NaN or infinite
            return xabs * yabs;
        } else if (Math.getExponent(max) - Math.getExponent(min) > EXP_DIFF_THRESHOLD_2D) {
            // value is completely dominated by max; just return max
            return max;
        }

        // compute the scale and rescale values
        final double scale;
        final double rescale;
        if (max > LARGE_THRESH) {
            scale = SCALE_DOWN;
            rescale = SCALE_UP;
        } else if (max < SAFE_SCALE_UP_THRESH) {
            scale = SCALE_UP;
            rescale = SCALE_DOWN;
        } else {
            scale = 1d;
            rescale = 1d;
        }

        // initialise sum and compensation using scaled x
        final double sx = xabs * scale;
        double sum = sx * sx;
        double comp = DD.twoSquareLow(sx, sum);

        // add scaled y
        final double sy = yabs * scale;
        final double py = sy * sy;
        comp += DD.twoSquareLow(sy, py);
        final double sumPy = sum + py;
        comp += DD.twoSumLow(sum, py, sumPy);
        sum = sumPy;

        return Math.sqrt(sum + comp) * rescale;
    }

    /** Computes the Euclidean norm.
     * This implementation handles possible overflow or underflow.
     *
     * @param x first input
     * @param y second input
     * @param z third input
     * @return \(\sqrt{x^2 + y^2 + z^2}\)
     *
     * @see #L2
     * @see #EUCLIDEAN
     * @see #of(double,double,double)
     */
    private static double euclidean(final double x,
                                    final double y,
                                    final double z) {
        final double xabs = Math.abs(x);
        final double yabs = Math.abs(y);
        final double zabs = Math.abs(z);

        final double max = Math.max(Math.max(xabs, yabs), zabs);

        // if the max is not finite, then one of the inputs must not have
        // been finite
        if (!Double.isFinite(max)) {
            // let the standard multiply operation determine whether to
            // return NaN or infinite
            return xabs * yabs * zabs;
        }

        // compute the scale and rescale values
        final double scale;
        final double rescale;
        if (max > LARGE_THRESH) {
            scale = SCALE_DOWN;
            rescale = SCALE_UP;
        } else if (max < SAFE_SCALE_UP_THRESH) {
            scale = SCALE_UP;
            rescale = SCALE_DOWN;
        } else {
            scale = 1d;
            rescale = 1d;
        }


        // initialise sum and compensation using scaled x
        final double sx = xabs * scale;
        double sum = sx * sx;
        double comp = DD.twoSquareLow(sx, sum);

        // add scaled y
        final double sy = yabs * scale;
        final double py = sy * sy;
        comp += DD.twoSquareLow(sy, py);
        final double sumPy = sum + py;
        comp += DD.twoSumLow(sum, py, sumPy);
        sum = sumPy;

        // add scaled z
        final double sz = zabs * scale;
        final double pz = sz * sz;
        comp += DD.twoSquareLow(sz, pz);
        final double sumPz = sum + pz;
        comp += DD.twoSumLow(sum, pz, sumPz);
        sum = sumPz;

        return Math.sqrt(sum + comp) * rescale;
    }

    /** Computes the Euclidean norm.
     * This implementation handles possible overflow or underflow.
     *
     * @param v input values
     * @return \(\sqrt{v_0^2 + ... + v_{n-1}^2}\).
     *
     * @see #L2
     * @see #EUCLIDEAN
     * @see #of(double[])
     */
    private static double euclidean(final double[] v) {
        // sum of big, normal and small numbers
        double s1 = 0;
        double s2 = 0;
        double s3 = 0;

        // sum compensation values
        double c1 = 0;
        double c2 = 0;
        double c3 = 0;

        for (int i = 0; i < v.length; ++i) {
            final double x = Math.abs(v[i]);
            if (!Double.isFinite(x)) {
                // not finite; determine whether to return NaN or positive infinity
                return euclideanNormSpecial(v, i);
            } else if (x > LARGE_THRESH) {
                // scale down
                final double sx = x * SCALE_DOWN;

                // compute the product and product compensation
                final double p = sx * sx;
                final double cp = DD.twoSquareLow(sx, p);

                // compute the running sum and sum compensation
                final double s = s1 + p;
                final double cs = DD.twoSumLow(s1, p, s);

                // update running totals
                c1 += cp + cs;
                s1 = s;
            } else if (x < SMALL_THRESH) {
                // scale up
                final double sx = x * SCALE_UP;

                // compute the product and product compensation
                final double p = sx * sx;
                final double cp = DD.twoSquareLow(sx, p);

                // compute the running sum and sum compensation
                final double s = s3 + p;
                final double cs = DD.twoSumLow(s3, p, s);

                // update running totals
                c3 += cp + cs;
                s3 = s;
            } else {
                // no scaling
                // compute the product and product compensation
                final double p = x * x;
                final double cp = DD.twoSquareLow(x, p);

                // compute the running sum and sum compensation
                final double s = s2 + p;
                final double cs = DD.twoSumLow(s2, p, s);

                // update running totals
                c2 += cp + cs;
                s2 = s;
            }
        }

        // The highest sum is the significant component. Add the next significant.
        // Note that the "x * SCALE_DOWN * SCALE_DOWN" expressions must be executed
        // in the order given. If the two scale factors are multiplied together first,
        // they will underflow to zero.
        if (s1 != 0) {
            // add s1, s2, c1, c2
            final double s2Adj = s2 * SCALE_DOWN * SCALE_DOWN;
            final double sum = s1 + s2Adj;
            final double comp = DD.twoSumLow(s1, s2Adj, sum) +
                c1 + (c2 * SCALE_DOWN * SCALE_DOWN);
            return Math.sqrt(sum + comp) * SCALE_UP;
        } else if (s2 != 0) {
            // add s2, s3, c2, c3
            final double s3Adj = s3 * SCALE_DOWN * SCALE_DOWN;
            final double sum = s2 + s3Adj;
            final double comp = DD.twoSumLow(s2, s3Adj, sum) +
                c2 + (c3 * SCALE_DOWN * SCALE_DOWN);
            return Math.sqrt(sum + comp);
        }
        // add s3, c3
        return Math.sqrt(s3 + c3) * SCALE_DOWN;
    }

    /** Special cases of non-finite input.
     *
     * @param v input vector
     * @param start index to start examining the input vector from
     * @return Euclidean norm special value
     */
    private static double euclideanNormSpecial(final double[] v,
                                               final int start) {
        for (int i = start; i < v.length; ++i) {
            if (Double.isNaN(v[i])) {
                return Double.NaN;
            }
        }
        return Double.POSITIVE_INFINITY;
    }

    /** Computes the maximum norm.
     *
     * @param x first input
     * @param y second input
     * @return \(\max{(|x|, |y|)}\).
     *
     * @see #LINF
     * @see #MAXIMUM
     * @see #of(double,double)
     */
    private static double maximum(final double x,
                                  final double y) {
        return Math.max(Math.abs(x), Math.abs(y));
    }

    /** Computes the maximum norm.
     *
     * @param x first input
     * @param y second input
     * @param z third input
     * @return \(\max{(|x|, |y|, |z|)}\).
     *
     * @see #LINF
     * @see #MAXIMUM
     * @see #of(double,double,double)
     */
    private static double maximum(final double x,
                                  final double y,
                                  final double z) {
        return Math.max(Math.abs(x),
                        Math.max(Math.abs(y),
                                 Math.abs(z)));
    }

    /** Computes the maximum norm.
     *
     * @param v input values
     * @return \(\max{(|v_0|, \ldots, |v_{n-1}|)}\)
     *
     * @see #LINF
     * @see #MAXIMUM
     * @see #of(double[])
     */
    private static double maximum(final double[] v) {
        double max = 0d;
        for (final double d : v) {
            max = Math.max(max, Math.abs(d));
        }
        return max;
    }

    /**
     * @param a Array.
     * @throws IllegalArgumentException for zero-size array.
     */
    private static void ensureNonEmpty(double[] a) {
        if (a.length == 0) {
            throw new IllegalArgumentException("Empty array");
        }
    }
}